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The  spatio-temporal  variability  of  soil  and  climate  properties  make  the  pre- 
dictions of  solute  fate  and  transport  in  subsurface  formations  a  very  difficult  and  an 
uncertain  task.  The  recent  development  of  stochastic  analysis  in  the  field  of  hydrology 
has  been  mainly  focused  on  estimating  uncertainty  of  flow  and  transport  predictions 
under  steady  flow  conditions  in  the  saturated  zone  due  to  the  spatial  variability  of 
the  formation  properties.  This  research  attempts  to  quantify  the  prediction  uncer- 
tainty of  flow  and  transport  in  heterogeneous  soils  subject  to  transient  water  flow 
regime.  The  objectives  here  are  to  derive  and  test  a  one-dimensional,  stochastic, 
advective-dispersive  solute  transport  model  that  incorporates  the  variability  of  the 
rainfall,  saturated  hydraulic  conductivity,  linear  equilibrium  adsorption  coefficient, 
and  water  uptake. 
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First,  a  stochastic,  unsaturated  flow  model  is  developed  to  predict  approximate 
expressions  for  the  means  and  autocovariances  of  the  soil-water  content,  Darcy's  flux 
and  pore-water  velocity  as  a  function  of  the  temporally  random  rainfall,  spatially 
random  saturated  hydraulic  conductivity,  and  spatially  variable,  temporally  random 
water  uptake.  Then,  these  statistics  are  tested  against  extensive  synthetic  data  gen- 
erated from  Monte-Carlo  simulations  that  are  assumed  to  represent  field  measure- 
ments. Second,  using  the  pore-water  velocity  results,  a  one-dimensional,  stochastic, 
advective-dispersive  solute  transport  model  is  derived  in  a  Lagrangian  framework  to 
predict  approximate  closed-form  expressions  for  the  mean  and  variance  of  the  solute 
flux  for  inert  and  reactive  solutes  affected  by  retardation  resulting  from  a  spatially 
random,  linear  equilibrium  adsorption  process.  Then,  these  transport  statistics  are 
tested  against  Monte-Carlo  experiments  for  two  soils  types  and  rainfall  patterns.  Both 
the  numerical  Monte-Carlo  experiments,  and  the  one-dimensional  stochastic,  unsat- 
urated water  flow  model  demonstrate  the  dominance  of  the  rainfall  variability  over 
that  of  the  saturated  hydraulic  conductivity  in  determining  the  overall  uncertainty 
in  the  vadose  zone  for  the  soil  and  climate  combinations  studied  here.  In  compari- 
son to  the  Monte-Carlo  simulations,  the  predicted  flow  statistics  are  of  good  quality, 
particularly  when  the  variance  of  the  rainfall  is  small  and  the  soil  has  a  fine  texture. 
The  stochastic  solute  transport  model  also  yields  reasonable  estimates  of  the  mean 
solute  flux  breakthrough  curves  (BTC's)  but  tends  to  underestimate  the  standard 
deviation  of  the  BTC's  at  the  mean  center  of  mass.  The  additional  effect  of  water 
uptake  on  solute  transport  is  not  only  to  reduce  the  mean  and  the  variance  of  the 
solute  velocity,  but  also  to  increase  significantly  the  temporal  correlation  scale  of  the 
solute  velocity.  Therefore,  the  incorporation  of  water  uptake  results  in  smaller  and 
delayed  peak,  and  larger  spreading  of  the  mean  and  standard  deviation  of  the  solute 
flux  BTC's.  Compared  to  an  inert  solute,  consideration  of  the  adsorption  coefficient 
leads  to  smaller  mean,  but  larger  variance,  of  the  apparent  solute  velocity,  resulting 
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in  a  smaller  and  delayed  peak  and  larger  spreading  of  the  solute  flux  BTC's.  Practical 
implications  and  applications  of  these  theoretical  findings  are  discussed. 
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CHAPTER  1 
BACKGROUND 

1.1  Introduction 

During  the  next  century,  modern  agriculture  will  have  to  reconcile  simultane- 
ously two  opposite  challenges:  providing  food  for  a  growing  world  population  while 
maintaining  the  quality  of  the  ecosystems.  This  difficult  achievement  will  require  ef- 
ficient management  of  water  resources.  Because  modern  agriculture  relies  extensively 
on  water  use  and  chemical  application,  the  issue  of  maintaining  a  safe  and  plentiful 
food  and  water  supply  becomes  necessarily  an  issue  of  water  resource  management. 

The  philosophy  of  water  resources  management  must  be  comprehensive.  Be- 
cause water  is  a  limited  and  vulnerable  resource  shared  by  multiple  shareholders,  this 
philosophy  is  expected  to  include  economic,  social  and  scientific  aspects.  The  ulti- 
mate goal  of  water  resources  management  is  to  find  an  optimal  solution  of  water  usage 
given  a  set  of  constraints  established  for  each  aspect.  Given  the  complexity  of  this 
holistic  approach,  our  interest  here  focuses  only  on  the  scientific  aspect.  This  latter 
aspect  deals  with  optimizing  the  quantity  and  quality  of  water  given  some  physical 
constraints.  This  implies  the  need  to  conceive,  test  and  develop  techniques  to  better 
maintain  and  protect  the  quality  as  well  as  to  secure  the  quantity  of  water. 

Groundwater  is  a  valuable  water  resource  while  at  the  same  time  quite  vul- 
nerable. For  example,  more  than  half  of  the  population  (52.5%)  of  the  United  States 
relies  upon  a  groundwater  source  (Solley  et  al,  1988).  Furthermore,  where  the  sur- 
face water  supply  becomes  insufficient  or  is  of  deteriorating  quality,  industrial  as  well 
as  agricultural  activities  exercice  increasing  pressure  on  groundwater  because  of  its 
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high  quality  and  abundance.  Consequently,  threats  of  overpumping  and  of  contam- 
ination sources  impose  stresses  on  a  hydrologic  system  which  is  well  characterized 
by  its  slow  replenishing  rate  and  long  time  of  recovery.  As  a  result,  subsurface  hy- 
drologists  face  the  challenge  of  developing  forecasting  management  tools  in  order  to 
predict  the  quality  as  well  as  the  quantity  of  groundwater. 

The  challenge  of  modeling  groundwater  flow  and  solute  transport  leads  to  at 
least  two  difficulties:  i)  accurate  estimation  of  recharge  and  contaminant  discharge 
coming  from  the  vadose  zone  and  entering  the  saturated  zone,  and  ii)  detailed  de- 
scription in  space  and  time  of  the  hydraulic  and  chemical  properties  of  underground 
formations  at  the  scale  of  interest.  Recently  developed  stochastic  models  of  flow  and 
transport  have  tackled  issues  related  to  the  second  difficulty  and  improved  predictions 
compared  to  the  deterministic  models.  However,  as  a  consequence  of  the  high  level  of 
complexity  of  the  vadose  zone  flow  and  transport  mechanisms,  the  first  difficulty  has 
received  less  attention  than  the  second.  Therefore,  the  objective  of  this  current  disser- 
tation is  to  develop  and  test  an  innovative  stochastic  model  of  unsaturated  transient 
water  flow  and  inert  and  reactive  solute  transport  in  the  vadose  zone. 

1.2    Deterministic  Modeling  of  Flow  and  Transport 

From  the  simplest  to  the  most  complex,  there  are  three  categories  of  deter- 
ministic models  commonly  used  to  simulate  the  flow  and  transport  in  underground 
formations.  The  increasing  level  of  complexity  implies  that  better  predictions  can  be 
obtained  in  compensation  for  an  increased  requirement  for  larger  collection  of  input 
data. 

Conceptual  models  are  those  that  retain  some  of  the  physical  laws  (e.g.,  con- 
servation of  mass)  in  their  mathematical  formulation,  without  trying  to  be  exact 
representations  of  the  physical  reality.  They  are  commonly  based  on  analogies  of 
river  channels  and  basins  as  a  set  of  storage  reservoirs  with  varied  properties.  A  good 
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example  is  the  model  developed  by  Nash  (1957)  that  represents  a  river  basin  as  a 
set  of  linear  reservoirs.  This  type  of  model  is  frequently  used  to  simulate  large  scale 
hydrological  processes. 

The  second  category  includes  black-box  models  that  do  not  attempt  to  describe 
the  detailed  physical  processes.  These  models  are  based  on  empirical  relationships 
derived  from  stimulus-response  data.  The  black-box  is  considered  as  a  filter  through 
which  the  input  data  are  converted  into  outputs.  Such  models  are  well  suited  to 
treat  systems  that  are  linear  and  invariant  and  about  which  a  limited  amount  of 
information  is  available.  The  outputs  are  often  described  as  a  convolution  integral 
of  the  input  function  multiplied  by  an  impulse  response  function  that  is  defined  as 
the  filter  response  to  an  instantaneous  pulse.  The  validation  of  such  models  requires 
large  amount  of  input  and  output  data  so  that  the  impulse  response  function  can  be 
identified. 

The  third  category  includes  mechanistic  physical  models.  Since  these  models 
are  the  most  accurate  and  also  widely  used,  it  is  worthwhile  to  describe  them  in  de- 
tail. These  models  are  based  on  the  assumption  that  the  underground  formations  can 
be  represented  by  a  pseudo-continuum  with  a  characteristic  Representative  Elemen- 
tary Volume  (REV)  over  which  partial  differential  equations  of  flow  and  transport 
are  applied.  The  three-dimensional  mechanistic  flow  equation  of  water  in  variably 
saturated  conditions  is  established  in  term  of  the  pressure  head,  ip,  as  (Freeze  and 
Cherry,  1979)  : 

(CW  +  5s)^  =  V-[/C(^Vh]  (1.1) 

where  C{ip)  is  the  slope  of  the  soil  moisture  characteristic  release  curve  [L_1] 
defined  as  JJ,  Ss  denotes  the  specific  storage  [Zr1]  that  is  a  function  of  the  porous 
medium  compaction,  and  of  the  fluid  compressibility,  }C(ip)  is  the  hydraulic  conductiv- 
ity [LT-1]  usually  represented  by  a  symmetric  tensor  of  second  order,  tp  is  the  pressure 


4 


head  [L]  and  h  is  the  total  head  [L]  defined  as  ip+z.  Assuming  an  incompressible  fluid 
and  nondeformable  porous  medium  away  from  the  water  table  boundary,  the  unsat- 
urated flow  equation  can  be  rewritten  in  terms  of  the  volumetric  soil-water  content, 
9,  as  the  independent  variable  (Bhattacharya  et  al,  1976): 


where  V{6)  is  the  soil- water  diffusivity  [L2T_1]  defined  as  K{0)$.  The  deter- 
ministic approach  assumes  that  i)  each  hydraulic  property  of  the  governing  equation 
is  defined  at  the  scale  of  the  representative  elementary  volume  which  is  larger  than 
the  pore-scale  heterogeneity  (assumed  to  be  between  10"5  and  10"1  m),  and  is  known 
in  a  continuous  way  over  space  and  time,  ii)  the  initial  and  boundary  conditions  are 
also  perfectly  known  over  space  and  time,  and  iii)  the  flow  regime  follows  the  linear 
Darcy's  Law  which  is  valid  for  a  Reynolds  number  between  1  and  10  (Bear,  1972). 

The  classical  three  dimensional  transport  equation  or  so  called  Advection- 
Dispersion  Equation  (ADE)  is  established  as  (Sposito  et  al,  1986): 


—  =  V  ■  [BD(6)VC]  -  qVC  (1.3) 
dt 

where  C  is  the  resident  concentration  in  the  aqueous  phase  [ML-3],  q  is  the 
Darcy's  flux  vector  [LT~l]  defined  in  terms  of  9  as  -2?(0)V0  +  K(6),  and  D  is  the 
dispersion  tensor  representing  the  hydrodynamic  dispersion  coefficients  [L2T~1].  In 
addition  to  the  assumptions  stated  above  for  the  mechanistic  flow  model,  the  mech- 
anistic transport  model  assumes  that  the  hydrodynamic  dispersion  is  similar  to  a 
diffusive  process  and  thereby,  is  mathematically  represented  by  the  first  Fickian  law 
of  diffusion. 

Given  particular  boundary  and  initial  conditions,  these  mechanistic  models  of 
flow  and  transport  can  be  solved  analytically  or  numerically.  However,  except  under 


some  very  simplified  systems  where  analytical  solutions  can  be  obtained,  numerical 
analysis  is  usually  performed  to  solve  these  equations.  This  observation  is  especially 
true  for  unsaturated  flow  where  the  hydraulic  conductivity  is  a  non-linear  function 
of  the  dependent  variable  (pressure  head  or  soil-water  content)  and  thus  obtaining 
analytical  solutions  is  almost  impossible.  The  usual  alternative  is  to  approximate 
the  solution  using  the  finite  differences  or  finite  elements  methods.  Although  both 
numerical  approaches  are  conceptually  different,  they  ultimately  lead  to  a  set  of  N 
equations  with  N  unknowns  which  are  solved  directly  or  iteratively. 

Despite  their  differences,  these  three  categories  of  models  hold  the  following 
common  and  major  assumptions.  All  the  hydraulic  properties  of  the  underground 
formations  or/and  all  the  initial  and  boundary  conditions  have  to  be  described  with 
certainty  and  in  a  continuous  way  over  space  and  time.  It  is  important  to  emphasize 
that  these  assumptions  are  fundamental  for  the  mechanistic  models.  Such  modeling 
expectations  are  often  satisfied  in  laboratory  experiments,  but  are  not  relevant  for 
most  hydrological  applications. 

Observations  of  geologic  formations  indicate  i)  a  high  level  of  variability,  and 
ii)  a  structured  arrangement  of  the  hydraulic  properties.  This  is  particularly  true 
for  the  saturated  hydraulic  conductivity  (Freeze,  1975;  Gelhar  and  Axness,  1983). 
Furthermore,  intensive  collections  of  field  data  over  large  areas  are  not  an  alternative 
for  economical  and  computational  reasons.  Therefore,  in  addition  to  the  modeling 
error  due  to  the  mathematical  formulation  of  the  physical  process,  one  can  anticipate 
i)  the  difficulty  of  extrapolating  laboratory  results  to  field  scale  problems  character- 
ized by  high  spatio-temporal  heterogeneities,  and  ii)  that  the  quality  and  quantity  of 
information  required  to  estimate  input  parameters  is  crucial  to  making  reliable  predic- 
tions. The  limitations  of  the  deterministic  modeling  approach  have  given  incentives 
to  develop  an  alternative  method:  the  stochastic  modeling  approach. 
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1.3    Stochastic  Modeling  of  Flow  and  Transport 

When  the  available  collection  of  data  is  insufficient  to  characterize  all  the  vari- 
ability of  the  investigated  area,  the  stochastic  approach  is  justified  by  introducing 
a  probabilistic  framework.  Thus,  the  spatio-temporal  input  data  are  assumed  to  be 
represented  by  random  processes  or  random  fields.  In  contrast  to  a  single  random 
variable  which  is  described  completely  by  its  probability  density  function,  a  continu- 
ous random  field  requires  the  joint  probability  density  functions  between  all  points  in 
space  or  time  to  describe  the  process  completely.  Because  a  limited  amount  of  data 
is  usually  available,  it  is  virtually  impossible  to  determine  these  joint  densities  func- 
tions. Rather,  a  more  practical  approach  is  to  consider  only  the  first  two  ensemble 
moments  of  the  joint  probability  density  functions.  The  mean  or  first  moment  is  used 
to  represent  the  general  trend  of  the  data,  the  variance  is  used  to  represent  its  level 
of  variability,  and  the  covariance  structure  its  level  of  correlation  over  space  and/or 
time. 

In  order  to  estimate  ensemble  moments  from  single  realizations  of  soil  or  geo- 
logic properties,  stochastic  approach  requires  two  fundamental  assumptions:  station- 
ary and  ergodicity.  The  stationarity  assumption  states  that  the  ensemble  moments 
of  the  random  processes  are  invariant  in  space  and  time.  The  tools  of  geostatistics 
allow  assessment  of  the  validity  of  this  assumption  from  limited  field  measurements. 
The  ergodicity  assumption  postulates  that  the  sample  moments  estimated  from  a  sin- 
gle realization  approach  the  population  or  ensemble  moments  as  many  observations 
are  taken  over  a  large  domain.  This  second  assumption  implies  that  the  ensemble 
moments  are  stationary  and  that  the  domain  of  study  covers  multiple  correlations 
scale  to  insure  many  independent,  identically  distributed  samples. 

Assuming  that  the  input  data  are  random  processes  implies  that  the  outputs 
or  the  predictions  of  any  model  are  also  random.  For  example,  the  treatment  of 
the  saturated  hydraulic  conductivity  as  a  random  space  function  means  that  the 
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resulting  head,  Darcy's  flux,  pore-water  velocity  and  solute  concentration  are  also 
spatiotemporally  random  processes.  The  drawback  of  this  stochastic  method  lies 
in  the  difficulty  of  validating  the  stationary  and  ergodic  assumptions  for  the  input 
parameters.  A  review  of  current  stochastic  methods  are  presented  below. 

Monte-Carlo  simulations 

The  Monte-Carlo  technique  is  a  systematic  method  for  determining  the  statis- 
tical moments  of  a  population  using  statistical  inference.  The  moments  of  the  outputs 
are  calculated  from  repeated  numerical  simulations  that  involve,  in  sequence,  i)  gen- 
erating random  fields  of  input  data,  ii)  solving  the  flow  and  transport  equations,  and 
iii)  evaluating  the  sample  statistics  of  the  output  data.  This  technique  is  powerful 
and  is  widely  used  because  its  implementation  is  easy  and  it  can  address  complicated 
real  case  scenarios.  However,  its  application  is  restricted  by  its  high  computational 
cost  in  situations  where  numerical  stability  and  convergence  are  uncertain,  and  where 
statistical  inference  demands  a  large  number  of  realizations.  The  first  applications 
of  the  Monte-Carlo  simulation  in  the  field  of  hydrology  were  conducted  by  Freeze 
(1975),  Smith  and  Freeze  (1979a),  and  Smith  and  Freeze  (1979b)  for  saturated  flow 
conditions,  and  Unlu  et  al.  (1990)  and  by  Polmann  et  al.  (1991)  for  unsaturated  flow 
conditions. 

In  the  absence  of  detailed  data  from  large  field  experiments,  it  is  a  common 
practice  to  test  and  verify  the  ensemble  moments  of  flow  and  transport  derived  from 
analytical  developments  against  Monte-Carlo  simulations.  This  verification  has  been 
frequently  undertaken  in  the  saturated  zone  (Graham  and  McLaughlin,  1989a,  1989b). 
However,  investigations  in  the  vadose  zone  tend  to  use  Monte-Carlo  simulations  more 
to  understand  the  complexity  of  flow  and  transport  processes  (Russo,  1991,  1993a, 
1993b,  Russo  et  al,  1994a,  1994b,  1998)  rather  than  to  check  the  very  limited  number 
of  analytical  solutions. 
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Eulerian  perturbation  technique 

The  perturbation  technique  consists  of  decomposing  all  random  process  in  the 
governing  flow  and  transport  equations  as  sums  of  means  and  perturbations  around 
these  means.  As  a  result,  the  flow  and  transport  equations  are  decomposed  into  a 
mean  partial  differential  equation  and  a  series  of  coupled  stochastic  partial  differential 
equations.  Analytical  methods  and  numerical  techniques  have  been  used  to  solve  these 
equations.  The  spectral  analysis  is  an  analytical  method  that  has  been  widely  used  by 
making  the  following  assumptions,  i)  the  perturbation  is  a  stationary  random  process 
uniquely  represented  in  terms  of  the  Fourier-Stieltjes  integral  transform,  ii)  time  and 
space  dimensions  are  infinite,  and  iii)  the  variance  of  the  perturbation  is  small.  Gelhar 
(1974)  initiated  this  spectral  analysis  that  has  been  widely  applied  in  the  saturated 
zone  under  steady  state  flow  regime  for  stationary  process  (Bakr  et  al.  1978;  Gutjahr 
et  al,  1978;  Gelhar  and  Axness,  1983)  and  under  transient  flow  regimes  (Rehfeldt  and 
Gelhar,  1992).  The  spectral  technique  was  later  relaxed  using  the  Fourier  Transform 
for  nonstationary  process  (Li  and  McLaughlin,  1991,  1995),  and  under  transient  flow 
regimes  (Li  and  Graham,  1998).  Although  the  non-linear  behavior  of  the  hydraulic 
conductivity  makes  the  flow  and  transport  in  the  vadose  zone  much  more  complex,  a 
few  attempts  have  been  performed  using  the  spectral  analysis  under  steady-state  flow 
conditions  (Yeh  et  al.,  1985a,  1985b,  1985c)  and  also  under  transient  flow  conditions 
(Mantoglou  and  Gelhar,  1987a,  1987b,  1987c).  Numerical  techniques  have  also  been 
implemented  to  solve  the  coupled  stochastic  partial  differential  equations  for  flow 
and  tranport  under  saturated  steady  state  flow  conditions  ( Graham  and  McLaughlin, 
1989a,  1989b;  James  et  al.,  1998;  Li  and  Graham,  1998)  and  under  unsaturated 
steady-state  flow  (Zhang  et  al,  1998). 
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Lagrangian  perturbation  technique 

The  Lagrangian  perturbation  technique  was  initiated  by  Dagan  (1982,  1984) 
to  calculate  ensemble  moments  of  the  resident  concentration  in  saturated  heteroge- 
neous flow  conditions.  In  this  formalism,  the  solute  is  viewed  as  a  set  of  infinites- 
imal particles  flowing  through  a  set  of  steamtubes  representing  flow  paths  through 
the  heterogeneous  porous  medium.  A  particle  initially  located  at  the  vector  posi- 
tion x0  at  time  t0,  moves  in  a  steamtube  along  a  Lagrangian  trajectory  coordinate, 
x  =  Xt(t;  t0,  x0)  which  can  be  decomposed  in  terms  of  a  mean  and  a  perturbation: 

Xt(t,  x0,  t0)  =  X(t\  x0,  h)  +  Xd(t,  to)  (1-4) 

where  Xd  stands  for  the  displacement  perturbation  due  to  pore  scale  disper- 
sion. The  particle  trajectory  is  related  to  the  Eulerian  velocity  through  this  funda- 
mental kinematic  differential  equation: 

^  =  Vt  =  Ut(Xt,  t)  +  u(Xt,  t)  +  vd  (1.5) 
at 

where  Vt  is  the  Eulerian  particle  velocity  decomposed  in  term  of  a  mean,  Ut, 
a  perturbation  u  and  vd  a  Brownian  motion  velocity  associated  with  the  pore  scale 
dispersion.  Once  Xt  is  known,  the  resident  concentration  in  a  streamtube  can  be 
expressed  by: 

dC(x,t,xQ,t0)  =  —5{x-Xt{t;x0,t0))  (1.6) 
n 

where  dC  stands  for  a  differential  increment  of  concentration,  n  represents  the 
porosity,  dM  =  C0nd  x0  denotes  the  initial  mass  of  an  infinitesimal  solute  particle  in 
a  volume  dV0  =  ndx0  with  a  concentration  of  C0,  and  5  is  the  Dirac  function.  The 
concentration  at  a  position  x  and  time  t  over  a  volume  V0  becomes: 
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C(x,  t)=       <J[x  -  Xt(t,  x0,  t0)]Co(x0,  *o)  dx 


0 


(1.7) 


where  C0(x0,  t0)  is  the  initial  solute  concentration.  Assuming  a  normal  proba- 
bility density  function  for  the  random  particle  displacement  Xt,  Dagan  (1984)  derived 
approximate  expressions  for  the  particle  displacement  variance,  and  the  mean  and 
variance  of  the  resident  concentration  in  a  steady-state  flow  under  a  uniform  mean 
hydraulic  head  gradient. 

Dagan's  work  was  later  extended  to  account  for  reactive  solute  (Bellin  et  al., 
1993;  Bellin  and  Rinaldo,  1995),  to  examine  the  effect  of  recharge  on  flow  nonunifor- 
mity  and  macrodispersion  (Rubin  and  Bellin,  1990),  to  predict  solute  flux  (Cvetkovic 
et  al.  1992;  Dagan  et  al,  1992),  and  to  examine  the  effect  of  transient  flow  regimes 
on  transport  statistics  (Bellin  et  al,  1996;  Zhang  and  Neuman,  1996).  This  approach 
offers  the  advantage  of  deriving  analytical  solutions  for  transport  under  a  steady-state 
flow  regime  that  provide  excellent  quantitative  tools  for  hydrologists.  Despite  the  ele- 
gance of  the  mathematical  formulation,  this  method  does  requires  some  computational 
efforts  for  the  numerical  integration  of  the  second  order  ensemble  moments. 

Transfer  function 

The  development  of  a  stochastic  transport  model  using  the  transfer  function 
formalism  was  initiated  by  Jury  (1982)  and  then  further  generalized  by  Jury  et  al. 
(1986).  The  rate,  Qout(t),  at  which  a  solute  is  leaving  a  soil  transport  volume  is 
established  by  a  convolution  integral: 


Jo 

where  Qin  is  the  rate  at  which  the  solute  mass  is  entering  the  transport  vol- 
ume, and  g(Tufe\Tin)  is  the  solute  lifetime  or  travel  time  probability  density  function 


(1.8) 
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describing  the  distribution  of  solute  lifetimes  THfe  (or  travel  time)  conditional  on  the 
time  Ti„  at  which  solute  enters  the  transport  volume.  The  lifetime  distribution  repre- 
sents the  net  effect  of  all  the  physical,  chemical  and  biological  processes  taking  place 
during  the  solute  travel  within  the  transport  volume.  It  is  conditioned  on  the  input 
time  in  order  to  account  for  the  time-dependency  of  any  physical,  chemical,  and  bi- 
ological changes  that  may  affect  the  lifetime.  This  model  has  been  simplified  for  the 
steady-state  case  by  removing  the  conditioning  on  rin  and  developed  for  reactive  and 
degrading  solute  (Jury  and  Roth,  1990).  Its  application  under  transient  flow  regime 
has  been  successfully  tested  against  some  Monte-Carlo  simulations  by  Jury  and  Gru- 
ber  (1989)  where  the  cumulative  amount  of  net  applied  water  and  the  soil  variability 
were  random,  and  against  some  field  measurements  by  Jury  et  al.  (1990)  where  the 
cumulative  water  input  was  random. 

The  assumptions  and  limitations  of  transfer  function  models  have  been  dis- 
cussed in  an  earlier  section  dealing  with  black-box  models.  However,  it  is  important 
to  recall  that  the  quality  of  the  impulse  response  function  is  crucial  for  reliable  pre- 
dictions and  that  it  is  usually  difficult  to  extrapolate  results  from  one  site  to  another. 

Solute  flux  theory 

An  alternative  representation  of  the  solute  transport  problem  is  to  characterize 
the  solute  flux  concentration  distribution  rather  than  the  local  resident  concentration 
distribution.  In  one  dimension,  the  resident  concentration  is  defined  as  the  mass  of 
solute  in  the  aqueous  phase  per  unit  volume,  whereas  the  solute  flux  concentration 
is  defined  as  the  ratio  of  the  solute  mass  flux  to  the  water  flux  past  a  particular 
boundary.  Compared  to  the  resident  concentration,  the  solute  flux  is  a  more  real- 
istic measurement  since  it  does  not  assume  that  the  solute  particles  are  immobile. 
Therefore,  over  the  past  years,  determination  of  solute  flux  ensemble  statistics  has 
received  considerable  attention.  The  primary  works  of  Shapiro  and  Cvetkovic  (1988) 
and  Cvetkovic  and  Shapiro  (1990)  initiated  a  stochastic  analysis  of  the  solute  flux 
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theory  that  was  later  generalized  by  Dagan  et  al.  (1992)  and  Cvetkovic  et  al.  (1992). 
This  theory  is  based  on  a  Lagrangian  formulation  of  solute  transport. 

The  solute  is  represented  by  discrete  and  indivisible  particles  that  move  along 
independent  streamtubes  that  represent  the  flow  lines.  For  an  instantaneous  injection, 
the  mass  flow  rate,  S(t),  passing  across  a  control  plane  at  distance  x  from  the  starting 
position  a  is  formalized  by: 

S(t,x)  =  /    C0(a)6(t-r)da  (1.9) 
Jo 

where  r  is  a  random  travel  time  reflecting  the  uncertainty  associated  with  the 
variability  of  the  hydraulic  properties,  a  denotes  the  coordinates  of  a  point  within  V0. 
C0 (a)  is  the  spatially  varying  concentration,  and  6  denotes  the  Dirac  function.  Then, 
taking  the  first  and  the  second  centralized  moment  of  the  previous  equation  after 
assuming  an  one-particle  probability  density  function  and  a  two-particle  probability 
density  function  of  the  travel  time  leads  to  the  mean  and  the  variance  of  S(t,x), 
respectively.  The  statistics  of  the  travel  time  are  derived  and  calculated  from  those 
of  the  particle  displacement  or  from  those  of  the  Eulerian  velocity. 

It  should  be  noted  that  the  solute  flux  theory  and  the  transfer  function  model 
are  very  similar.  The  transfer  function  method  uses  impulse  response  functions  cal- 
ibrated to  input-output  data  while  the  solute  flux  theory  usually  tries  to  relate  the 
impulse  response  function  to  physical  processes  and  properties  of  the  soil.  The  so- 
lute flux  theory  has  been  used  to  treat  transport  of  reactive  solute  in  saturated  flow 
conditions  (Dagan  and  Cvetkovic,  1996;  Cvetkovic  and  Dagan,  1996),  transport  in 
the  vadose  zone  (Destouni,  1992a,  1992b;  Destouni,  1993),  and  also  transport  in  the 
soil-groundwater  system  (Destouni  and  Graham,  1995). 

The  solute  flux  theory  does  not  face  the  scale  effect  of  the  deterministic  ADE 
equation.  It  provides  elegant  analytical  solutions  that  are  useful  tools  for  hydrologists. 


However,  it  does  require  specifying  the  single  as  well  as  the  joint  probability  density 
function  of  the  travel  time,  which  are  not  easily  determined. 

1.4    Outline  of  the  Dissertation 

As  stated  earlier  in  the  introduction,  the  current  work  focuses  on  predict- 
ing the  timing  and  distribution  of  water  and  solute  flux  through  the  heterogeneous 
vadose  zone  subject  to  a  temporally  random  flow  regime.  So  far,  most  previous  in- 
vestigations dealing  with  the  flow  and  transport  statistics  have  assumed  steady-state 
flow  (Dagan  and  Bresler,  1979;  Yeh  et  al,  1985a,  1985b,  1985c;  Destouni,  1992a, 
1992b;  Destouni,  1993;  Destouni  and  Graham,  1995),  or  have  run  a  limited  num- 
ber of  Monte-Carlo  simulations  under  transient  flow  regime  (Russo,  1991;  Russo  et 
al,  1994a,  1994b,  1998).  However,  because  the  Monte-Carlo  simulations  are  com- 
putationally costly,  and  little  research  has  been  reported  on  exact  or  approximate 
analytical  solutions  for  flow  and  transport  statistics  in  unsaturated  conditions  sub- 
ject to  transient  flow  regime  (Mantoglou  and  Gelhar,  1981a;  Jury  et  al.  1990a),  there 
is  a  need  to  develop  and  validate  new  analytical  methods  for  this  purpose. 

This  dissertation  will  develop  a  stochastic  analysis  of  flow  and  solute  transport 
in  heterogeneous  soils  subject  to  temporally  random  rainfall  and  spatially  variable, 
temporally  random  water  uptake.  Both  inert  and  sorbing  solutes  will  be  consid- 
ered. Chapter  2  contains  an  evaluation  of  the  effect  of  random  rainfall  and  random 
saturated  conductivity  on  the  statistics  of  the  soil- water  content,  Darcy's  flux  and 
pore-water  velocity  in  a  one-dimensional  vertical  soil.  Semi-analytical  solutions  of 
the  flow  statistics  are  derived  from  the  one-dimensional  Richards'  equation.  Intensive 
Monte-Carlo  simulations  for  two  types  of  soils  and  rainfall  patterns  are  performed  to 
test  the  validity  of  the  stochastic  unsaturated  flow  model. 

Chapter  3  extends  the  study  of  Chapter  2  by  incorporating  a  spatially  distributed 
temporally  random  water  uptake  into  the  governing  flow  equation,  and  the  effect  of 
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water  uptake  on  flow  statistics  is  examined.  Chapter  4  analyzes  the  solute  flux  statis- 
tics for  an  inert  and  a  reactive  solute  under  flow  regimes  developed  in  Chapters  2 
and  3.  Closed  form  analytical  solutions  of  the  solute  flux  statistics  are  derived  in  the 
solute  flux  framework  that  relies  on  Lagrangian  transport  theory.  The  effects  of  the 
variabilities  of  the  soil,  rainfall,  water  uptake,  and  linear  equilibrium  sorption  on  the 
mean  and  standard  deviation  of  the  solute  flux  are  examined.  Intensive  Monte-Carlo 
simulations  for  two  types  of  soils  and  rainfall  patterns  are  performed  to  test  the  va- 
lidity of  the  stochastic  advective-dispersive  transport  model  for  inert  solute  transport 
under  a  transient  flow  regime. 


CHAPTER  2 

UNSATURATED  TRANSIENT  STOCHASTIC  FLOW  MODEL 

2.1  Introduction 

Accurate  descriptions  of  contaminant  transport  in  soils  are  essential  for  man- 
agement of  hazardous  wastes.  Therefore,  reliable  predictions  of  solute  transport 
through  soils  requires  understanding  and  accurately  modeling  the  flow  mechanisms 
taking  place  in  transient,  variably  saturated,  porous  media.  This  task  is  hindered  for 
three  reasons:  i)  the  nonlinearity  of  the  flow  equation,  ii)  the  spatial  variability  of  the 
soil  properties,  and  iii)  the  spatio-temporal  variability  of  the  weather  patterns. 

The  application  of  stochastic  methods  integrates  the  spatio-temporal  variabil- 
ity of  the  hydrological  processes  in  a  probabilistic  modeling  framework.  Except  for 
some  very  simplified  cases,  such  as  the  assumption  of  a  steady-state  flow  regime, 
analytical  solutions  to  multi-dimensional,  stochastic,  unsaturated  flow  and  transport 
problems  are  not  possible.  Thus,  numerical  solutions  are  the  rule  rather  than  the 
exception.  Although  the  steady-state  assumption  is  not  always  realistic  in  the  va- 
dose  zone,  this  simplification  has  been  widely  used  because  it  offers  the  advantage  of 
examining  only  the  effect  of  soil  spatial  variability  on  predicting  flow  and  transport 
statistics  under  unsaturated  flow  conditions  (Dagan  and  Bresler,  1979;  Yeh  et  ai, 
1985a  ,  1985b  ,  1985c  ;  Destouni,  1992a;  Russo,  1993a  ,  1993b). 

Single,  large-scale  realization  or  numerical  Monte-Carlo  simulations  have  been 
frequently  used  to  examine  transient  unsaturated  zone  processes  (Russo  et  al.,  1994a, 
1994b,  1998).  However,  few  attempts  at  deriving  analytical  expressions  for  unsatu- 
rated flow  ensemble  moments  under  transient  regime  have  been  reported  so  far.  The 
most  comprehensive  works  are  those  of  Mantoglou  and  Gelhar  (1987a,  1987b,  1987c). 
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Mantoglou  and  Gelhar  developed  a  framework  for  modeling  large-scale  transient  un- 
saturated flow  in  spatially  variable  soil.  They  used  the  three-dimensional  Richards 
equation  to  obtain  a  partial  differential  equation  for  the  mean  and  perturbation  of 
pressure  head.  The  resulting  equations  were  solved  using  spectral  analysis  in  order  to 
predict  mean  pressure  head  and  pressure  head  variance  in  stratified  soils.  Polmann  et 
al.(1991)  tested  the  theoretical  developments  of  Mantoglou  and  Gelhar  (1987a,  1987b, 
1987c)  by  running  a  single  numerical  experiment. 

Parlange  et  al.  (1992)  formulated  a  first-order  autoregressive  Markovian  model 
for  the  soil-water  content  on  the  basis  of  a  lumped  hydrologic  budget  where  the 
applied  water  events  were  white  noise  in  structure.  They  successfully  tested  the 
statistics  of  the  soil- water  content  under  uniform  soil  conditions  against  neutron  probe 
measurements.  Katul  et  al.  (1997)  used  a  semi-analytical  method  to  derive  the 
mean  and  the  time-space  variance  of  the  soil-water  content  by  accounting  for  the 
spatial  interaction  between  soil-water  content  and  two  other  processes,  root  uptake 
and  vertical  flux  gradient.  Their  results  successfully  validated  the  mean  soil- water 
content  against  some  laboratory  measurements. 

These  previous  investigations  show  limitations  in  the  validation  of  their  theo- 
retical development  and  in  their  possible  application.  The  works  of  Mantoglou  and 
Gelhar  (1987a,  1987b,  1987c)  use  spectral  analysis  which  requires  local  stationarity 
of  random  processes,  and  they  do  not  account  explicitly  for  any  temporally- varying 
term  forcing  the  system  to  be  transient.  Thus,  their  framework  makes  their  findings 
difficult  to  be  applied  to  real  life  scenarios  that  are  usually  bounded  in  time  and  space. 
Their  work  needs  to  be  further  tested  since  only  one  numerical  experiment  has  been 
performed  over  a  short  and  unique  rainfall  event.  The  model  developed  by  Parlange 
et  al.  (1992)  is  a  lumped  model  and  therefore  does  not  allow  prediction  of  soil-water 
content  over  depth.  The  work  of  Katul  et  al.  (1997)  was  validated  against  laboratory 
measurements  and  numerical  experiments  that  considered  the  saturated  hydraulic 
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conductivity,  the  porosity  and  the  near-saturated  pressure  head  as  random  processes. 
However,  their  analysis  did  not  include  temporally  varying  boundary  conditions.  All 
these  aformentioned  works  focus  on  predicting  statistics  of  soil- water  content;  they 
do  not  address  predictions  of  other  important  flow  parameters  such  Darcy's  flux  and 
pore-water  velocity. 

From  the  foregoing  review,  it  is  apparent  that  compared  to  the  spatial  vari- 
ability of  soil  properties,  the  incorporation  of  temporally  varying  surface  infiltration 
processes  in  predicting  vadose-zone  flow  ensemble  moments  has  received  little  at- 
tention in  either  numerical  simulation  or  analytical  development.  Therefore,  in  this 
Chapter  the  statistics  of  the  soil-water  content,  Darcy's  flux  and  pore-water  velocity 
will  be  derived  which  explicitly  account  for  temporal  variability  of  the  rainfall  pattern 
and  spatial  variability  of  the  saturated  hydraulic  conductivity.  The  objectives  of  this 
Chapter  are  i)  to  develop  a  one-dimensional  stochastic  model  of  unsaturated  flow  in 
heterogeneous  soil  subject  to  a  transient  flow  regime,  ii)  to  examine  the  relative  effects 
of  the  rainfall  uncertainty  versus  soil  variability  on  the  ensemble  moments  of  soil  wa- 
ter content,  Darcy's  flux  and  pore-water  velocity,  and  iii)  to  test  the  developed  model 
against  some  quasi-realistic,  two-dimensional,  numerical  Monte-Carlo  experiments. 

2.2    Analytical  Development 
2.2.1    Non-dimensional  Form  of  the  Richards'  Equation 

In  unsaturated  soils  subject  to  rainfall  and  irrigation,  water  flows  mainly  in 
the  vertical  direction.  This  assumption  is  supported  by  field  experiments  (Butters 
et  al,  1989)  and  numerical  simulations  (Russo,  1991).  Taking  advantage  of  this 
simplification,  a  one-dimensional  equation  for  transient  unsaturated  flow  is  assumed 
in  which  the  saturated  hydraulic  conductivity  and  the  rainfall  have  random  structure 
over  space  and  time,  respectively: 


18 


-  =  ¥z{Kt{z)V{B)  -)  ^— 


(2.1) 


where  Ka{z)V{9)  is  the  soil-moisture  diffusivity,  K{9)  is  the  dimensionless  rel- 
ative hydraulic  conductivity  which  is  a  function  of  the  local  volumetric  soil-water 
content  9,  and  Kt{z)  is  a  spatially  random  saturated  hydraulic  conductivity.  Note 
that  in  equation  (2.1),  the  Darcy's  flux  oriented  positively  downward  along  the  ver- 
tical direction  z  is  defined  as: 


q  =  -Kt(z)V(9)^-+Ka(z)K{B)  (2.2) 
oz 

In  order  to  account  for  soil  variability,  a  lognormal  distribution  is  assumed  for 
the  saturated  hydraulic  conductivity,  Ks*{z)  =  K*Gefr**)  where  KG*  is  the  geometric 
mean  [LT~l],  /  is  a  constant  parameter  equivalent  to  the  standard  deviation  of  In 
/C*  and  r\  is  a  normally  distributed  random  space  function  of  zero  mean  and  unit 
variance.  It  is  further  assumed  that  all  these  quantities  are  relevant  at  the  scale  of 
the  representative  elementary  volume,  and  that  the  air  phase  is  immobile.  In  this 
analysis,  equations  (2.1)  and  (2.2)  are  made  nondimensional  by  using  the  following 
scaling  factors: 


z  =  —,t=  — V  =  — ,  Ka(z)  =  -£V-,  and?  =  — 
hb-  hi  hi  K-g  *"G 


where  hi  is  the  local  bubbling  pressure  [L]  and  the  asterisks  indicate  the  di- 
mensional quantities.  Neglecting  surface  runoff  and  water  uptake  due  to  evapotran- 
spiration,  and  assuming  a  deep  water  table,  the  dimensionless  initial  and  boundary 
conditions  are  established  for  a  semi-infinite  domain  as: 
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r)9 

-fCs(z)V(9)  —  +  )Cs(z))C(0) 


71  +  eS(t) 


(2.3) 


z=0,t 


89 


-Ks(z)V{9)  —  +  lCs(z)K{9) 


2=00,  t 

d\z,t=0 


n  (2.4) 

Oi  (2.5) 

where  TZ  is  the  dimensionless  mean  rainfall  defined  as         (q0*)  denotes  the 
expected  value  of  the  rainfall  series,  e  is  a  small  dimensionless  parameter  defined  as 
with  <79-  representing  the  standard  deviation  of  the  rainfall  series,  8(t)  is  a  random 
function  of  time  with  a  zero  mean  and  unit  standard  deviation  and  Oi  is  the  initial 
soil-water  content.  The  ()  stands  for  the  expected  value. 

2.2.2    Power  Series  Solution 

The  difficulty  of  solving  the  partial  differential  equation  for  unsaturated  flow 
rests  in  the  non-linear,  random  coefficients  and  boundary  conditions.  Here,  we  assume 
that  the  solution  to  equation  (2.1)  can  be  approximated  by  a  power  series  solution  of 
the  form  {Mei,  1997): 


9(z,t;e,f)  =  J2J2d^z^elf3  t2"6) 

i=0  j=0 

where  small  values  of  e  and  /  insure  the  convergence  of  the  series.  An  arbitrary 
non-linear  function  g(9)  can  be  expanded  about  #0o  in  a  Taylor's  series  as  follows: 


9(0)  =  g(9oo)  + 


dg 
09 


(9  -  900)  + 


0oo 


dO2 


(9  -  900f 


+ 


Boo 


where  each  of  the  partial  derivatives  is  to  be  evaluated  at  900.  Note  that  each 
term  (0  —  #oo)n  can  be  simplified  as: 


(0  -  900)n  =  ^^(MK/^-fl 

,i=0  j=0 


00 


,i=l  j=l 
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The  exponential  form  of  the  spatially  random  log-saturated  hydraulic  conduc- 
tivity can  be  approximated  in  a  Taylor  series  expansion  as: 


f2  2 

e}r}  ~  1  +  frj  +  tJj-  +  .  • 


Substituting  these  series  expansion  in  equations  (2.1)-(2.5),  the  left  hand  side 
in  equation  (2.1)  becomes  to  second  order: 


89  _  d{900  +  eQl0  +  f  ggi  +  ef  9n  +  e2  O20  +  P  O02  +  ■  •  ■ ) 
dt  ~  dt 

and  the  right  hand  side  of  equation  (2.1)  becomes: 


(2.7) 


d_ 

dz 


/V 


(    %)+SL  (£^io  +  /^oi  +  €/eu+  \ 


d0  \0W 

e2e20  +  f2e02  +  ...) 


,  1  <Pv 
~t~2  ae2 


(eOw  +  f  Ooi  +efOn  + 12  920+ 

V  pe02  +  . ..)'  +  ...  J 

(1  +  frj+  §-z  (0OO  +  e9l0  +  f  901  +  e/  »u  +  <?  #20+ 

+  ...  /2^02  +  ---) 

/    £(0oo)+  §L  (e*u>  +  /*tt+«Mi+  \ 
e202o  +  /20o2  +  -.-) 
{ee10  +  fe0i  +ef6n  +  e2  920+ 

V  /20o2  +  -.-)2  +  ---  7 
Substituting  the  power  series  expansion  of  0  in  equations  (2.3)-(2.5),  the 

boundary  conditions  become 


■  1  d±K 

"1"2  ae2 


Ooo 
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+  . . . 


/^(M+f  L  (^lO  +  ^Ol+X 
efOu  +  e2  92Q  +  f2  902  +  ...) 


+ 


1  92"D 


2  ae2 


\ 


(€010  +  /  #01  + 

ef8n+e2e20+ 

f2e02  +  ...)2  + 


J 


£{000 +  €010  +  / 001+ C/ 011+ 

e202o  +  /20O2  +  .-O 
/  /C(0oo)+  f  L00  (€010  +  / 0oi+  \ 
e/0n+e2  02O  +  /20O2  +  ---) 

(e^io  +  /0oi  + 


,i  a^c 
"^2  ae2 


<?oo 

e/  0ii  +  e2  020+ 
V  /20o2  +  -..)2  +  ---  /  J 


+  ... 


/ 


2>(*oo)+§ L  (^io  +  /  0oi+ 
e/0n  +  e2  02o  +  /20o2  +  -.-) 

(6^1O  +  /0Ol+€/011  + 


a2^ 
2  ae2 


7e  +  e5(2,t) 
at  2  =  0 


V       e2  920  +  f2  602  +  ...)2  +  ...  J 

&  (000  +  €  0io  +  /  001  + 

e/0n+e2  02O  +  /20O2  +  ---) 
/    /C(0oo)+f  L00  (e0io  +  /0oi+  \ 
ef6n+e2  620  +  f2  0O2  +  --O 
+J|f      (e0lo  +  /0oi+e/0n+ 

^00 

V       e202O  +  /20O2  +  ...)2  +  ---  / 


while  the  initial  condition  is 


71 

at  z  =  oo 


(000  +  €  010  +  /  001  +  €/  011  +  e2  02O  +  /2  002  +  •  ■  • )  L>t=0    =  0i 


The  power  series  solution  of  0(z,  t)  requires  solving  the  governing  flow  equation 
and  its  initial  and  boundary  conditions  for  each  up  to  i+j  equal  to  the  desired 
order  of  approximation  (which  is  typically  not  higher  than  second  order).  Collecting 
terms  of  O(e°f0)  leads  to  the  governing  equation  of  the  zero  order  soil-water  content, 
0oo : 
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^O0-I:[^oo)^f +  Wo)]  =  o 


dt  dz' 


d800 


(-mo)^+W) 

an 

^oo|2it=0  =  0» 


71 


z=0,i 


z=oo,t 


n 


(2.8) 
(2.9) 

(2.10) 
(2.11) 


Note  that  this  is  the  non-linear  Richards  equation  which  must  in  general  be 
solved  numerically.  However,  the  boundary  conditions  considered  here  are  constant  in 
space  and  time.  Therefore,  away  from  the  initial  condition,  the  steady-state  solution 
to  this  equation  requires  that: 


00, 


-2>(0Oo)^f+  K{eQ0)  =  n 


(2.12) 


at  all  locations  in  the  soil  profile.  A  valid  solution  for  this  equation  is  obtained 
by  setting  the  zero  order  soil-water  content  gradient  to  be  zero,  ^  =  0,  so  that 
IC(900)  =  TZ.  This  equation  can  then  be  solved  for  600  once  a  relative  hydraulic 
conductivity  function  has  been  specified. 

Collecting  terms  of  0(ef°),  and  noting  that  ^f-  =  0  leads  to  the  governing 
equation  for  9W: 


dt 


dz"2 


d)C 
~d0 


86 


10 


(-v(eoo)—+  — 

(    r>(()  \ddl°  +  dlC 


Boo 

0io ) 


dz 


=  0 


Boo 


2  =  0,t 


2  =  +  CX),t 


e 


iolz,t=o 


0 


(2.13) 
(2.14) 

(2.15) 
(2.16) 
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This  equation  resembles  an  advection-dispersion  equation  with  constant  co- 
efficients (obtained  by  evaluating  the  appropriate  dimensionless  soil-water  content 
characteristic  relationships  at  the  solution  to  0Oo)  and  a  zero  mean,  unit  variance 
temporally  correlated  random  boundary  forcing  term,  5(t).  Since  the  expected  value 
of  the  forcing  term  is  zero,  the  expected  value  of  6W  will  also  be  zero. 

Collecting  terms  of  0(e°f),  rearranging  the  equation  so  the  lower  order  forcing 
terms  appear  on  the  right-hand  side,  and  noting  that  ^  =  0  leads  to  the  governing 
equation  for  0Oi: 


d9oi     V(R  /0oi  dJC 


000 


d90l 
dz 


0 


d6Qi  ,  dJC 
86 


+ 


Boo 


001  ) 


z=0,t 


(-^00)-^-+  -88eoo6oi) 


2  =  0O,t 

d0l\z,t=0 


dz 
-r]JC(9oo) 

-r]K,(900) 
0 


(2.17) 
(2.18) 

(2.19) 
(2.20) 


This  equation  resembles  an  advection-dispersion  equation  with  constant  co- 
efficients and  both  a  random  interior  and  random  boundary  forcing  term  driven  by 
T](z).  Since  the  expected  value  of  all  forcing  terms  are  zero,  the  expected  value  of  this 
equation  will  also  be  zero. 

Collecting  terms  of  0(el fl)  and  rearranging  the  equation  so  that  the  lower 
order  forcing  terms  appear  on  the  right-hand  side,  and  noting  that  =  0  leads  to 
the  governing  equation  for  On. 
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, .  .d29n  dK 

at    v{9oo)~d^ +  ~w 


Ooo 


dOn 

dz 


d_ 

dz 


dV\      a  d6io 

ae  \e00  °01  dz 
av\     n  in 
ae  1 0oo  10  9z 
a2K 


ae2 


#01010  + 


(-  Wo)1f+  f  L>  ) 


z=0,t 


000 

Woo)1f- 

/   dV\      a    Mlll-L  \ 

wle00yio  9,  + 

#01  #10  + 

r/Woo)1f+ 

dK  I       /)  1 
dV       n  dtfoj. 

ae  1 0oo  10 


d2K 
d92 


z—oo,t 


d2K 

ae2 


Ooo 


#01#10  + 

,30, 


V 


dJC\ 

ae  \e0, 


1  >i  y 


=  0 


(2.21) 


(2.22) 


(2.23) 


(2.24) 


#n|2,t=o 

This  equation  also  resembles  an  advection-dispersion  equation  with  constant 
coefficients  and  both  a  random  interior  and  random  boundary  forcing  term  driven  by 
901, 910,  and  r](z).  Note  that  the  expected  value  of  this  equation  is  different  from  zero, 
because  the  expected  values  of  its  forcing  terms  are  not  necessarily  zero. 

Collecting  terms  of  0(e2f°),  rearranging  the  equation  so  that  the  lower  order 
forcing  terms  appear  on  the  right-hand  side,  noting  that        =  0,  and  also  making 

io  bit 


use  of  the  fact  that  #10^  =  ^^,  leads  to  the  governing  equation  for  02O 


89 
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dt 


V{9m) 


d29 


20  dfC 
dz2  39 


89 


20 


(  V(9  )d02O+  dK 
(-V(900)—+  — 


#20  ) 


(-2>(M 


d9- 


20 


dz 


dK 
89 


Ooo 
#20  ) 


dz 


z=0,t 


z=+co,t 


8_ 

dz 

dV 
~89 

8V_ 
~89 

0 


av\  JLiia 
ae  1 0oq  dz  2 


a2fc 
ae2 


0oo 


62 
gio 

2 


d  92 


000 


_10  ,  SMC 
dz  2  d92 


d  92 


dz  2  d92 


000 


000 


(2.25) 


f  (2.26) 


f  (2.27) 
(2.28) 
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This  is  an  advection-dispersion  equation  with  constant  coefficients,  a  random 
interior  forcing  term,  and  a  random  boundary  forcing  term.  Note  that  the  expected 
value  of  this  equation  is  not  zero,  but  is  driven  by  (9f0)  =  ajl0. 

Collecting  terms  of  0(e°f2)  and  rearranging  the  equation  so  that  the  lower 
order  forcing  terms  appear  on  the  right-hand  side  after  noting  that  =  0,  and  also 
making  use  of  the  fact  that  0Oi  ^  =  ;  "af  >  we  have  : 


d902     „fa  ,d2902  dIC 


39, 


02 


000 


dz 


{'v{doo)~dV+  do 


#02  ) 


Ooo 


z=0,t 


(-V(900)—  +  — 


9o2 ) 


Boo 


z—+oo,t 


02|z,i=0 


dv\  d  £ai 
ae  \e00  dz  2 


+ 


*01  + 


ai/c 
d      962  e00 

dz  v(V(9oo)d-t+ 
dK\     ggoy  | 
Iffoo  9z 
.  ^-/C(^oo) 


/  22I     A?k+  \ 

90  I0oo  9z  2 


ae2 


2  ~ 


0oo 

V[V(900)d-t+ 

V  £/C(0oo) 
/  ^2 1  A^oi_(- 


30  1 6»oo  3*  2 

a2^  fgi.  , 
a02  fl    2  ^ 


\ 


=  0 


(2.29) 


(2.30) 


(2.31) 


(2.32) 


This  is  an  advection-dispersion  equation  with  constant  coefficients,  a  random 
interior  forcing  term,  and  a  random  boundary  forcing  term.  Because  the  expected 
value  of  its  forcing  term  is  different  from  zero,  the  expected  value  of  this  equation 
is  also  not  zero.  However,  note  that  there  are  no  transient  forcing  terms  in  this 
equation.  Therefore,  away  from  the  initial  condition  a  temporally  stationary  steady- 
state  solution  will  exist. 
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2.2.3    Governing  Equation  of  the  Flow  Ensemble  Moments 
Soil-water  content  statistics 

The  soil-water  content  is  a  spatio-temporally  random  process  which  depends 
upon  the  random  forcing  terms  driving  the  flow  system.  In  the  absence  of  detailed 
measurements  of  these  forcing  terms,  it  is  impossible  to  characterize  9(z,  t).  However, 
the  central  tendency  of  its  behavior  can  be  estimated  by  its  mean  value  and  the 
expected  variability  around  that  mean  by  the  variance.  Taking  the  expected  value  of 
equation  (2.6)  leads  to  the  following  equation  for  the  mean  soil-water  content: 


(9(z,  t))  =  (900)  +  €  (9W)  +  e2  (0O2>  +  ef  (On)  +  f  (&oi)  +  f  (flm)  +  •  ■ 


(2.33) 


Equations  (2.8)  -  (2.11)  are  deterministic,  therefore  (900(z,t))  =  90o-  Fur- 
thermore, as  mentioned  previously,  since  the  expected  values  of  the  forcing  terms  of 
equations  (2.13)  and  (2.17)  are  zero,  (6W),  and  (0Oi)  are  also  zero.  Assuming  that  5 
and  r)  are  uncorrected  ((0Oi0io)  =  0),  the  expected  value  of  equation  (2.21)  leads  to 
(On)  =  0. 

Similarly,  taking  the  expected  value  of  equations  (2.25)  -  (2.27),  and  noting 
that  (920)  =  00lo  we  have: 


d(920)     n(a  ,d2{920)  6K 


d(920) 
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(-V(90Q) 


d(920)  ,  dlC 
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d2K 

dO2 
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dV  1 

ae  \ 

a  ^.n  | 
0OO  dz    2  1 

d2K 

ae2 

e00  2 

<*»>Uo   =  0 
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Thus,  (020 )  is  the  solution  to  a  constant  coefficient  advection  dispersion  equa- 
tion forced  both  within  the  domain  and  on  the  boundary  by  a  linear  operator  on  <x|M 
which  is  potentially  a  function  of  time  and  space. 

Similarly  taking  the  expected  value  of  equations  (2.29)  -  (2.31),  noting  that 


8Bqi(x)\  _ 


(#oi)  =  cr|oi ,and  again  noting  that  (viz)-^) 


9pr,em  (fV) 
dz 


we  have: 
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dz 
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Boo 
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90  Loo  ^W(2>  z) 
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^K(0oo) 


ffP  I 
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86  1 0oo  92   2    T  ae2 
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^00  J  ^  

96  Loo  P^oii*''2) 
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9oo 


z'=z 


+ 
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av\     JL^Li.  _i_  ^ 

90  1 0oo  9z    2    T  902 


£>(0oo) 


9Pr?9m  (z'.x) 


9z 


000 
z'=z  + 


90  Loo  ^>r,9°1       ^  'z'=z  + 


(0O2)L>t=o  =  o 

Thus,  (0O2 )  is  the  solution  to  a  constant  coefficient  advection  dispersion  equa- 
tion forced  both  within  the  domain  and  on  the  boundary  by  a  linear  operators  on 
a!  ,  and  Pve01(z',  z).  Since  away  from  the  initial  conditions,  the  forcing  terms  are  not 
a  function  of  time,  a  temporally  stationary  steady-state  solution  to  (0O2)  exists. 

Utilizing  these  developments,  the  second  order  mean  soil- water  content  can  be 
rewritten  as: 


(0(z,t))    =   0oo  +  e2(02oOM;0oo,<702lo))  + 

/2  (002  {Z,  t]  0oo,  00O1 ,  P 70oi 

(z',z)))  (2.34) 
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To  a  second  order  the  soil- water  content  autocovariance,  9,  Pg6,  can  be  esti- 
mated as  follows: 

Pee(z,t;z',t')  =  ((*(*,*)  -  (*(*,*)»  W  ft  ~  W  ft))) 

,17  9oo{z,t)  +  edw(z,t)  +  f0Oi (*, *)  +  •••-  \ 

~Hv  (Mm)  +  --.)  / 

"V  (M*',*')  +  ■••)  /J 

=  e2  (M*.  *)  M*'.  O)  +  f  (^oi  <)  0oi  (**,  f  )>  +  £/  (M*.  t)  601  (*',  t')) 
=  e2  P9lofllo  (*,  t;  z',  t1)  +  f2  Pg0l  g0l  (*>  t;  z\  t')  +  ef  Pg10  g01  (z,  t;  z',  if) 

(2.35) 

Assuming  no  correlation  between  the  rainfall  and  the  saturated  hydraulic  con- 
ductivity, the  third  term  on  the  right  hand  side  of  equation  (2.35)  is  zero.  Thus, 
the  approximation  of  the  soil-water  content  covariance  indicates  that  the  contribu- 
tion of  the  rainfall  and  soil  variability  are  additive.  The  components  of  this  sum  can 
be  evaluated  from  equations  (2.13)  -  (2.16)  and  (2.17)  -  (2.20)  in  a  variety  of  ways. 
The  option  chosen  here  is  to  develop  partial  differential  equations  for  Pel0s(z,  t\  z',  t'), 
Pg10gl0(z,t;z',t'),  PBoi11(z,t;z',t')  and  P6oigol  (z,  t;  z1,  f )  from  equations  (2.13)  -  (2.16) 
and  (2.17)  -  (2.20),  respectively. 

Pre-multiplying  equations  (2.13)-(2.16),  by  a  soil-water  content  perturbation, 
Oio{z',t'),  at  location  z'  and  time  t'  independent  of  z  and  t  and  taking  the  expected 
value,  leads  to  the  governing  equation  of  Pg10gl0 , 
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dPeme10(z',t'-z,t)  d2Pei0eJz',t';z,t)  dK 

 at  v[doo)       d?       +  ~m 


dPe10e10(z',t';z,t) 
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dz 


=  0 
(2.36) 
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dz  89 
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Pe109io\z',t';z,t=0  =  ® 


PgloS(z',t';z  =  0,t) 


z',t';z=0,t 


=  0 


z'  ,t';z—+oo,t 


Again  this  equation  resembles  an  advection-dispersion  equation  with  constant 
coefficients,  driven  by  the  boundary  forcing  term  PgloS.  The  boundary  forcing  term 
may  be  evaluated  by  post-multiplying  equations  (2.13)-(2.16)  by  S(z  -  0,t),  at  an 
independent  location  z=0  and  a  time  t,  and  taking  the  expected  value, 
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z'=+oo,t';z=0,t 


In  a  similar  manner  the  equation  for  PeQ1e01  is  obtained  by  pre-multiplying 
equation  (2.17)-(2.20)  by  90i(z',t'),  at  an  independent  location  z'  and  a  time  t',  and 
taking  the  expected  value, 
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The  forcing  term  for  these  equations  may  be  evaluated  by  post-multiplying 
equation  (2.17)-(2.20)  by  r}(z),  at  a  location  z  independent  of  z',  and  taking  the 
expected  value, 


dP9oiT)(z',t';z)  ^PeoiV(z',t';z)  dK 
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=  -X^(0oo)^iw(2  > z) 
=  -K.{900)Pm(z',z) 


z' =+ 00, t';z,t 


(2.43) 

(2.44) 
(2.45) 


Once  P6lo0lo,  Peoi0ol  and  Pfloi7?  have  been  solved  a210  a2Q1,  and Pe0lT1(z,  z)  can  be 
evaluated  and  the  second  order  corrections  to  the  mean  water  content,  (9u(z,  t;  0Oo, 
Pgoir)(z',z))),  (92o(z,t;9oo,<J2810{z,t))),  and  {9Q2(z,  t;  90Q,  aj01  (2),  PfloilJ(z',  2)))  can  be 
determined. 

Darcy's  flux  statistics 

The  spatio-temporal  variability  of  the  specific  discharge  or  Darcy's  flux,  from 
the  vadose  zone  to  the  saturated  zone  is  driven  by  spatially  variable  soil  properties  as 
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well  as  by  spatio-temporal  variability  of  surface  infiltration.  In  the  absence  of  detailed 
descriptions  of  soil  properties  and  rainfall  patterns,  it  is  impossible  to  characterize 
the  Darcy's  flux.  However,  its  general  behavior  may  be  summarized  using  its  first 
and  second  moments. 

Inserting  the  power  series  expansion  of  9,  and  the  Taylor's  series  expansion  of 
V(0),  JC(6)  and  efr>  into  equation  (2.2)  for  the  Darcy's  flux,  and  noting  that  =  0 
leads  to  the  following  second  order  equation: 
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Taking  the  expected  value  of  equation  (2.46)  leads  to  the  following  equations 
for  the  means  Darcy's  flux: 


(?)  =  /c(M 


Zero  and  first  order 


1  dV\ 


(<?)=/C(M+e2[-2  88\800  dz 
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1  dV 

2  86  \600  dz 


Second  order 


do2 


where  the  expressions  (flio^1)  and  (#oi^)  simplify  to  and  \ 


8*1 


dz  I  \^ui  dz  I  *^    ^  2    dz  2    dz  ' 

respectively,  the  component  —  T>(90Q)^-+  || |0oo  #oi+?7£(0oo)  equals  to  zero  (equation 
(2.18)  )  and  the  cross-covariance  between  r\  and  5,  P^,  is  assumed  to  be  null.  Note 
that  the  second  order  mean  Darcy's  flux  can  be  calculated  once  cr|io,  oj  ,  (0n),  (#20) > 
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and  Pvg0l(z  =  z')  are  known.  The  second  order  Darcy's  flux  autocovariance  can  be 
approximated: 

Pqq(z,t;z',t')  =((q(z,t)  -  (q(z,t)))  (q(z',t')  -  (q(z',t')))) 

=([(  Ww)  +  e[-V(6w)^ +^91^,1)]-^)  ) 
•  (  ^00)+e[-2>(doo)^^+  §  I^M^Ol-W  )]) 


2jT)lQ  ,2d2Peloel0(z',t';z,t)  dK 


)2^>fioeio(2''*';2'i) _ 


,  dIC      fdP0weMAz,t)  dPg10g10(z',t';z,t)^ 
W  e00  dz  dz  (2.47) 

where  again  it  was  noted  that  -P(0oo)^f  +  %\Boo  #oi  +  vWoo)  =  0  and  it 
was  assumed  that  PStl  =  0.  Note  that  due  to  the  one-dimensional  problem  and  the 
selection  of  the  boundary  conditions  the  autocovariance  of  the  Darcy's  flux  depends 
only  upon  the  variability  of  the  rainfall. 

Pore-water  velocity  statistics 

Infiltrating  water  carries  inert  and  reactive  solute  through  soils.  Quantify- 
ing the  level  of  solute  mass  entering  the  saturated  zone  requires  knowledge  of  the 
pore-water  velocity  in  the  vadose  zone.  Therefore,  the  uncertainty  of  solute  transport 
prediction  relies  on  estimating  the  spatio-temporal  variability  of  the  pore-water  ve- 
locity due  to  the  spatially  random  structure  of  the  saturated  hydraulic  conductivity 
and  the  temporally  random  structure  of  rainfall  series.  Our  interest  focuses  on  the 
mean  and  the  variance  of  the  pore-water  velocity,  v(z,t)  defined  as  the  ratio  of  two 
random  processes, 

»(M)  =  «  (2-48) 

In  order  to  derive  the  statistics  of  v,  a  series  expansion  of  equation  (2.48) 
about  #oo  is  required  (Mood  et  ai,  1974), 
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,    #v  <?(z,t)(0(*,t)-0oo)  ,  q{z,t){6(*,t)  -  6oq)2  ,  . 

V{z,t)  =  —  ^   +   Z3   +  •••  VA*> 

POO  pOO  ^00 

Collecting  term  of  order  e,  e2,  /,  /2,  and  ef  after  substituting  equations  (2.6) 
and  (2.46)  into  equation  (2.48)  leads  to 
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Evaluating  the  expected  value  of  velocity  at  the  zero,  first  and  second  order 
leads  to: 
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where  the  same  simplifications  and  assumptions  used  to  derive  the  expected 
values  of  Darcy's  flux  have  been  applied.  The  second  order  dimensionless  pore-water 
velocity  autocovariance  can  be  approximated: 
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(2.51) 


where  it  was  assumed  that  =  0.  Note  that  the  autocovariance  of  pore- 
water  velocity  depends  upon  both  the  variability  of  the  rainfall  as  well  as  that  of  the 
saturated  hydraulic  conductivity  due  to  the  dependence  of  the  variability  of  soil- water 
content  on  saturated  hydraulic  conductivity  variability. 
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2.2.4    Semi-analytical  Solutions  of  the  Flow  Ensemble  Moments 

Equations  (2.47)  and  (2.51)  indicate  that  the  ensemble  moments  Pqq  and  Pvv 
are  both  linear  functions  of  Pe10e10  and/or  Pe01e01-  Determining  Pe10e10  and  Pe0i60i 
requires  first  solving  for  PgloS  and  Peoin-  To  solve  the  equations  for  P6l0s,  a  dou- 
ble Laplace  transform  is  applied  over  t  and  t'  in  equation  (2.37)  and  its  respective 
boundary  and  initial  conditions, 
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where  the  double  tilde  represents  the  double  Laplace  transform  and  s  and  s' 
are  the  Laplace  variables  corresponding  to  t  and  t',  respectively.  The  solution  of  this 
second  order  linear  and  nonhomogeneous  differential  equation  is  given  by  Jury  and 
Roth,  (1990b,  p.  48) 


P6loS(z',  s';  z  =  0,s)=  —  —  Pss(s,  s')e^r{1-^  (2.52) 


with  &  «  Ji  + 

Similarly  performing  a  double  Laplace  transform  over  equation  (2.36)  for 


Pe10e10  and  its  associated  initial  and  boundary  conditions  leads  to: 
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The  solution  of  this  equation  is  similar  to  equation  (2.52)  and  equals  to: 


Pel0elQ{z\  s' ]  z  -  0,s)  -  gic,  . 


(2.53) 


with  f:  =    /l  +  4ff  °°|2 .  Substituting  equation  (2.52)  into  equation  (2.53) 
V      ( 98  kc/ 

gives  the  soil-water  content  autocovariance  function  in  the  Laplace  domain: 
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Equations  (2.42)  -(2.45)  for  P0oir?  and  equations  (2.38)  -(2.41)  for  P6oie01  have 
no  temporal  forcing  terms.  Therefore,  away  from  the  initial  condition,  the  governing 
equations  can  be  evaluated  at  steady-state: 
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and 
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The  forcing  terms  in  the  top  and  bottom  boundary  conditions  for  equations 
(2.55)  and  (2.56)  are  identical  and  also  identical  to  the  interior  forcing  term.  Thus,  if 
Pm  is  assumed  to  be  stationary,  an  infinite  domain  solution  can  be  found  that  satisfies 
these  equations.  Taking  the  double  Fourier  transform  of  equation  (2.55)  leads  to  an 
algebraic  expression  for  the  cross-spectrum  between  0Oi  and  77: 


HC(9Q0)STm(ki) 
-V{6m)kx  -  i  f  ^ 


(2.57) 


where  ki  is  the  wave  number  corresponding  to  the  spatial  lag,  f ,  of  the  sta- 
tionary function,  Pm.  Taking  the  double  Fourier  transform  of  equation  (2.56)  yields 
the  following  expression  for  the  autospectrum  of  60i. 
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Inserting  equation  (2.57)  into  equation  (2.58)  leads  to: 
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The  autocovariance,  Pe01e01  can  be  found  by  taking  the  inverse  Fourier  trans- 
form of  S0oleol: 
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Taking  the  double  Laplace  transform  of  equations  (2.35)  leads  to  the  soil-water 
content  autocovariance  in  the  Laplace  domain: 


P„(zf,  s';  z,  s)  =  e2P^0 (z',  s';  z,  s)  +  ^Pe01e01  M  *)  (2.60) 

where  Pe01e01(z';  z)  has  been  evaluated  in  equation  (2.54)  and  Pe01e0i  {z'\ z)  is 
time-invariant.  Taking  the  double  Laplace  transform  of  equations  (2.47)  leads  to  the 
Darcy's  flux  autocovariance  in  the  Laplace  domain  that  is  only  a  function  of  the 
temporal  component  of  the  soil- water  content,  Pe10ei0'- 

Pqq{z\  s';  z,  s)  =  e2[iW*',  s';  z,  s)'"^'      4  (2.61) 

Taking  the  double  Laplace  transform  of  equation  (2.51)  leads  to  the  pore- water 
velocity  autocovariance  in  the  Laplace  domain  which  depends  upon  both  components 
of  the  soil-water  content,  Pg10e10  and  Pe0le0l'- 
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where  it  is  noted  that  Pe01e01  is  time-invariant.  As  soon  as  the  autocovariances 
of  the  rainfall  and  saturated  hydraulic  conductivity  are  specified,  a  two-dimensional 
Fast  Fourier  Transform  algorithm  can  be  used  to  invert  equations  (2.60),  (2.61),  and 
(2.62)  to  find  P0g(z',f;z,t),  Pqq(z' ,f;  z,t),  and  Pvv(z' ,t';  z,t)  in  real  time.  Appendix 
1  summarizes  the  relationship  between  the  Laplace  and  Fourier  transforms. 
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2.3    Illustration  of  Theoretical  Results 

The  mean  and  the  autocovariance  of  the  soil-water  content,  Darcy's  flux  and 
pore-water  velocity  are  illustrated  for  a  heterogeneous  fine  sandy  soil  subject  to  a 
random  rainfall  pattern  characteristic  of  Tampa,  Central  Florida.  The  soil-moisture 
diffusivity  and  relative  hydraulic  conductivity  are  modeled  using  the  Brooks  and 
Corey  relationships  (Brooks  and  Corey,  1964).  This  empirical  model  is  well  adapted 
for  sandy  soils  (Lappala  et  al,  1993)  and  usually  fits  measured  field  data  better 
than  the  less  complex  Gardner- Russo  model  (Zhang  et  al,  1998).  The  Brooks-Corey 
model  represents  the  soil-moisture  diffusivity  and  relative  hydraulic  conductivity  as 
follows: 


vs  —  vr 

T>(9)  ^       ( ^ _  ^r  )2+^ 

X(9S  —  9r)  9S  —  9r 


where  9S  the  saturated  soil- water  content,  equivalent  to  the  porosity,  6r  the 
residual  soil-water  content,  and  A  the  grain  size  distribution  parameter.  The  Taylor's 
series  expansion  of  V  and  K,  about  0Oo  are  expressed  as: 
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The  soil-moisture  diffusivity  was  determined  using  the  first  term  of  the  Tay- 
lor's series  expansion  of  V.  Thus,  this  constant  value  corresponds  to  a  soil-moisture 
diffusivity  estimated  at  a  soil-water  content  averaged  over  space  and  time.  This  lin- 
earization differs  from  that  performed  by  Philip  (1969,  p.  260)  in  the  sense  it  does  not 
match  the  integral  behavior  of  the  infiltration  process  in  the  limit  as  t  ->  0  and  differs 
from  that  used  for  a  linear  soil  {Gardner,  1980)  in  the  sense  that  it  depends  upon 
the  soil- water  content.  The  quantity  ^  was  also  determined  at  a  time-and-space 
averaged  soil- water  content.  This  second  estimation  differs  from  that  used  by  Philip 
(1969,  p.  260)  in  the  sense  that  it  does  not  match  the  infiltration  behavior  as  t  -4  oo. 
Table  2.1  summarizes  the  deterministic  properties  of  the  soil  used  to  illustrate  the 
theoretical  results,  and  Table  2.2  summarizes  the  statistical  parameters  of  the  rainfall 
series  and  saturated  hydraulic  conductivity. 

The  variance  of  the  saturated  hydraulic  conductivity  was  selected  as  the  upper 
bound  of  /  so  that  it  would  insure  a  converging  power  series  solution  for  6.  The 
statistics  of  the  rainfall  were  calculated  from  long  term  weather  records  from  Tampa, 
Florida. 

2.3.1    Mean  of  the  Soil- water  Content,  Darcv's  Flux  and  Pore- water  Velocity 

After  specifying  the  relative  hydraulic  conductivity  relationship,  the  zero  order 
soil-water  content  is  determined  by  solving  equation  (2.12): 


^-(^Ff^+3  =  o  (2>63) 

v  s  '  "r 

This  nonlinear  equation  is  solved  using  the  Newton-Ralphson  method.  Then, 
the  first  order  dimensional  means  of  the  Darcy's  flux  and  of  the  pore-water  velocity 
are  easily  found  to  be  equal  to: 
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#00 

where  the  dimensionless  quantity  q  and  v  are  multiplied  by  the  scaling  factor 
K*G  to  obtain  the  dimensional  Darcy's  flux  and  pore-water  velocity.  Since  the  power 
series  of  6  converges  rapidly  due  to  the  small  values  of  e  =  1.9  x  10~5  and  /  =  1, 
the  calculation  of  the  second  order  mean  for  the  soil-water  content,  Darcy's  flux  and 
pore-water  velocity  has  not  been  undertaken. 

2.3.2    Autocovariances  of  the  Soil-water  Content.  Darcy's  Flux  and  Pore-water  Velocity 

Evaluating  the  autocovariances  of  the  soil-water  content,  Darcy's  flux  and 
pore-water  velocity  requires  specifying  a  temporal  correlation  structure  for  the  rain- 
fall, a  spatial  correlation  structure  for  the  saturated  hydraulic  conductivity,  and  ex- 
panding the  nonlinear  functions  of  soil-moisture  diffusivity  and  relative  hydraulic 

conductivity  around  900.   A  negative  exponential  function  is  assumed  for  the  sta- 

\t-t'\ 

tionary  temporal  covariance  of  the  rainfall  series,  PSs(t,t')  =  e  xs  where  Xs  is  the 
temporal  correlation  length  of  the  rainfall  series.  The  spatial  covariance  of  the  satu- 
rated hydraulic  conductivity  is  also  modeled  with  a  stationary  negative  exponential 
function,  Pw(£)  =  where  is  the  spatial  correlation  length  along  the  vertical 
direction  and  f  =  z  —  z'  the  spatial  lag. 

The  double  Laplace  transform  of  Pgs  is  given  by: 


PS5(z  =0,s;z  =  0,s)   =    (1  +  v/)(s  +  s0(1  +  M  (2-64) 

Substituting  equation  (2.64)  into  equation  (2.54)  yields  the  following  temporal 
component  of  the  soil-water  content  autocovariance  in  Laplace  space: 
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e-^rd-^)+-2^-(i-€l)  (2.65) 
The  Fourier  transform  of  Pm  is  given  by: 


where  fci  is  the  frequency  corresponding  to  the  spatial  lag  £.  Using  Mathe- 
matica  (Wolfram,  1991),  the  inversion  of  equation  (2.59)  leads  to: 


^oi»oi(0  =  \£(0oo)  (777   ax  I    \5     T>/fl  \2\  + 
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where  £  =  \z  -  z'\.  The  dimensional  autocovariances  of  the  Darcy's  flux  and 
pore-water  velocity  in  the  Laplace  domain  are  obtained  by  inserting  equation  (2.65) 
and  (2.67)  into  equations  (2. 61)  and  (2.62),  inverting  to  real  time  using  the  two- 
dimensional  Fast  Fourier  Transform  algorithm  and  then  multiplying  by  the  scaling 
factor  ICq2. 

2.3.3    Results  and  Discussion 

The  normalized  spatio-temporal  autocovariances  of  the  soil-water  content, 
Darcy's  flux,  and  pore-water  velocity  are  illustrated  in  Figures  2.1-2.12.  In  these 
figures  the  autocovariances  are  normalized  by  dividing  by  the  variance  calculated  at 
the  land  surface.  The  relative  effect  of  the  saturated  hydraulic  conductivity  and  the 
rainfall  variability  on  the  soil-water  content  autocovariance  is  illustrated  in  Figures 
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2.1  and  2.2.  These  figures  show  that  for  a  given  depth  the  soil- water  content  autoco- 
variance  does  not  decay  completely  to  zero  for  large  time  lags.  This  is  due  to  the  fact 
that  the  rainfall  driven  component  of  the  soil- water  content  autocovariance  decays  to 
zero  for  large  time  lags  while  the  hydraulic  conductivity  driven  component  is  evalu- 
ated at  a  zero  spatial  lag  that  yields  a  constant  value  equal  to  the  spatial  variance. 
Furthermore,  these  figures  indicate  that  at  shallow  depths,  the  effect  of  the  rainfall 
variability  dominates  over  that  of  the  saturated  hydraulic  conductivity.  However,  the 
importance  of  rainfall  variability  decreases  as  depth  increases. 

Figures  2.1  and  2.3  show  the  normalized  soil- water  content  autocovariance  at 
two  different  depths.  These  figures  indicate  that  away  from  the  initial  condition,  at 
a  given  depth  the  soil-water  content  variance  and  autocovariance  become  stationary 
with  time.  However,  the  soil-water  content  is  non-stationary  over  depth,  with  a  vari- 
ance decreasing  and  the  temporal  correlation  increasing  as  depth  below  land  surface 
increases.  Figures  2.1  and  2.3  show  also  that  the  plateaus  are  identical  indicating 
the  spatially  stationary  contribution  of  the  log-saturated  hydraulic  conductivity.  As 
expected,  these  figures  also  demonstrate  that  the  soil-water  content  autocovariance 
is  a  symmetric  function. 

The  temporal  behavior  of  the  normalized  soil- water  content  autocovariance  is 
explained  by  the  equations  governing  the  propagation  of  soil-water  content  through 
the  soil  profile.  The  linearization  of  the  0-based  Richards  equation  about  0Oo  trans- 
forms a  nonlinear  partial  differential  equation  of  flow  into  a  series  of  linear  advection- 
dispersion  equations.  Thus,  the  soil-moisture  diffusivity  acts  as  a  dispersive  term 
that  dampens  the  effect  of  the  random  rainfall  forcing  term  over  depth.  Near  the  soil 
surface,  the  soil-water  content  autocovariance  has  a  larger  variance  and  a  temporal 
correlation  length  of  the  same  magnitude  as  the  rainfall  correlation  length.  At  deeper 
depths  the  effect  of  soil-moisture  diffusivity  is  demonstrated  by  reducing  the  variance 
and  increasing  the  temporal  correlation  scale  of  the  soil-water  content  autocovariance. 
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The  deterministic  initial  condition  explains  the  time  delay  in  reaching  the  variance 
plateau  at  deeper  depths  (Fig.  2.3).  Figures  2.4  and  2.5  show  the  autocovariance 
of  the  Darcy's  flux  and  the  pore-water  velocity  at  a  depth  of  0.5  meter  normalized 
by  their  respective  variance  at  the  soil  surface.  The  behavior  of  these  two  autoco- 
variances  is  similar  to  that  of  Peg;  that  is  increasing  temporal  correlation  scale  and 
decreasing  variance  with  depth  below  land  surface.  It  is  important  to  notice  that 
the  contribution  of  the  soil  heterogeneity  is  zero  on  Darcy's  flux  autocovariance,  very 
small  on  the  pore-water  velocity  autocovariance  and  small  on  the  soil-water  content 
autocovariance,  as  indicated  by  the  plateau  value  in  Figures  2.1,  2.4,  and  2.5. 

Comparing  Figures  2.6  and  2.7  demonstrates  the  relative  effect  of  the  rainfall 
variability  and  of  the  soil  variability  on  the  spatial  evolution  of  the  soil- water  content 
autocovariance;  i.e,  the  soil-water  content  autocovariance  values  are  dominated  by 
the  rainfall  variance,  particularly  at  shallow  depths.  Note  the  stationarity  of  the  soil- 
water  content  autocovariance  due  to  only  the  log-saturated  hydraulic  conductivity, 
as  indicated  by  the  plateau  in  Figure  2.7.  Figures  2.6  and  2.8  show  the  normalized 
soil-water  content  autocovariance  over  depth  at  two  different  times.  These  figures 
confirm  that  the  soil-water  content  autocovariance  is  non-stationary  over  depth.  Like 
the  temporal  behavior  of  the  soil-water  content  autocovariance,  its  spatial  behavior 
is  also  controlled  by  the  flow  governing  equation.  The  stationary  rainfall  series  is 
responsible  for  the  identical  profiles  of  soil-water  content  autocovariance  at  different 
times.  As  depth  increases,  the  effect  of  the  soil-moisture  diffusivity  becomes  more 
evident,  the  variance  decreases  and  the  soil- water  contents  are  correlated  over  larger 
distances.  The  spatial  behavior  of  the  Darcy's  flux  and  of  the  pore-water  velocity 
autocovariance  exhibit  similar  behavior  as  that  of  the  soil-water  content  (Figures  2.9 
and  2.10).  Note  that  at  the  soil  surface  the  Darcy's  flux  variance  is  equivalent  to  the 
rainfall  variance. 
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Figures  2.11a,  b  and  2.12a,  b  describe,  in  time  and  space,  the  effect  of  the 
saturated  hydraulic  conductivity  and  rainfall  variability  {a2qQif  =  9.81  x  10_5ra2/d2) 
on  the  autocorrelation  of  the  soil-water  content  and  pore-water  velocity.  Figure  2.11a 
shows  that  i)  an  infinite  temporal  correlation  scale  of  the  soil  water  is  obtained  with 
no  rainfall  variability,  ii)  increasing  the  saturated  hydraulic  conductivity  variance 
increases  the  temporal  correlation  scale  of  the  soil- water  content  more  than  an  increase 
of  its  correlation  scale,  and  iii)  the  soil  variability  increases  the  temporal  correlation  of 
the  soil-water  content.  Therefore,  the  temporal  correlation  of  the  soil-water  content 
is  more  governed  by  the  variance  of  the  saturated  hydraulic  conductivity  than  by  its 
spatial  correlation  scale.  Figure  2.11b  shows  that  i)  the  largest  spatial  correlation  scale 
of  the  soil  water  is  obtained  with  no  rainfall  variability,  ii)  increasing  the  correlation 
scale  of  the  saturated  hydraulic  conductivity  increases  the  spatial  correlation  scale 
of  the  soil- water  content  more  than  an  increase  of  its  variance,  and  iii)  the  rainfall 
variability  decreases  the  spatial  correlation  of  the  soil- water  content.  Therefore,  the 
spatial  correlation  scale  of  the  soil- water  content  is  more  controlled  by  the  correlation 
scale  of  the  saturated  hydraulic  conductivity  than  by  its  variance.  The  effect  of  the 
soil  heterogeneity  on  the  pore-water  velocity  autocovariance  is  insignificant  (Figures 
2.12a,  b).  Because  of  the  form  of  the  bottom  boundary  condition,  the  statistics  of 
the  Darcy's  flux  are  not  affected  by  the  soil  variability  for  the  1-D  problem. 

The  main  findings  of  these  theoretical  results  are  summarized  as  follows.  Away 
from  the  initial  condition,  the  soil-water  content,  Darcy's  flux,  and  pore-water  veloc- 
ity become  asymptotically  stationary  in  time  because  of  the  temporal  stationarity  of 
the  rainfall  series.  However,  the  advection-dispersion-like  flow  equation  dampens  the 
effect  of  the  rainfall  series  over  depth,  making  the  soil-water  content,  Darcy's  flux 
and  pore-water  velocity  non-stationary  over  depth.  The  variability  of  the  soil-water 
content,  Darcy's  flux  and  pore- water  velocity  are  dominated  by  the  variability  of  the 
rainfall  at  shallow  depths.  The  importance  of  the  saturated  hydraulic  conductivity 
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variability  on  the  soil-water  content  variability  increases  away  from  the  land  surface, 
but  is  relatively  minor  for  the  pore-water  velocity  over  all  the  depths  studied.  It 
should  be  noted  that  the  solution  technique  used  here  does  not  require  a  continu- 
ous covariance  of  rainfall.  It  is  possible  to  take  the  Laplace  transform  of  a  discrete 
covariance  function  using  a  numerical  Fast  Fourier  Transform  algorithm. 

2.4    Validation  Against  Monte-Carlo  Simulations 

The  development  of  the  stochastic  unsaturated  flow  model  requires  many  sim- 
plifying assumptions  and  linearizations.  To  test  the  validity  of  these  assumptions, 
results  from  the  stochastic  unsaturated  flow  model  are  compared  against  intensive 
two-dimensional  Monte-Carlo  simulations  that  are  assumed  to  be  representative  of 
field  data.  Monte-Carlo  simulations  with  random  rain  but  uniform  soil  properties 
were  conducted  to  test  the  linearization  of  the  soil-moisture  diffusivity  and  relative 
hydraulic  conductivity  functions  around  the  zero  order  mean  value  of  the  soil-water 
content.  Monte-Carlo  simulations  with  both  random  rainfall  and  random  soil  prop- 
erties were  conducted  to  further  test  the  linearization  of  efr>  and  the  assumption  of 
one  dimensional  flow. 

2.4.1  Introduction 

The  numerical  experiments  consisted  of  a  Monte-Carlo  simulation  that  repeat- 
edly simulated  400  days  of  infiltration  of  rainfall  into  homogeneous  and  heterogeneous 
isotropic  unsaturated  soil.  The  domain  of  simulation  represents  the  vertical  cross  sec- 
tion of  a  30  m  wide  by  10  m  deep  soil  overlying  a  horizontal  water  table  (Fig.  2.13). 
Rainfall  enters  the  system  uniformly  over  the  entire  land  surface.  No  flow  occurs  on 
either  side  of  the  domain. 
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A  modified  version  of  the  variably  saturated  flow  and  transport  code,  VS2DT 
(Lappala  et  al,  1993,  Healy,  1990),  was  used  to  simulate  the  infiltration  experi- 
ment. This  code  solves  the  head-based  Richards  equation  for  flow  and  the  advection- 
dispersion  equation  for  transport.  The  Brooks  and  Corey  model  was  utilized  to  model 
the  soil-moisture  characteristic  curves.  The  code  was  modified  in  order  to  accept  ran- 
dom inputs  of  saturated  hydraulic  conductivity  and  random  time  series  of  rainfall 
fluxes  at  the  soil  surface.  Appropriate  spatial  and  temporal  discretization  were  se- 
lected to  insure  convergence,  stability  and  accuracy  of  the  numerical  scheme.  Vertical 
and  horizontal  grid  sizes  of  .25  m  and  a  maximum  time  step  of  .01  days  were  uti- 
lized to  satisfy  the  Peclet,  Neuman  and  the  Courant  number  conditions.  The  initial 
condition  was  established  using  the  Brooks  and  Corey  equilibrium  pressure  profile 
above  the  horizontal  water  table.  The  rainfall  was  treated  as  a  temporally  random 
boundary  conditions  at  the  land  surface. 

Intensive  Monte  Carlo  simulations  including  150  runs  each  were  performed  for 
two  rainfall  patterns  representative  of  Tampa,  (Florida)  and  Olympia  (Washington 
State)  and  two  soil  types  for  which  the  saturated  hydraulic  conductivity  was  assumed 
to  be  either  uniform  or  spatially  random.  The  random  fields  of  saturated  conductivity 
were  generated  using  the  turning  band  algorithm,  originally  developed  by  Journel  and 
Huijbregts,  (1978)  and  then  modified  by  Akpoji  (1993).  The  simulations  were  run 
for  a  fine  sandy  soil  (Brooks  and  Corey,  1964)  and  a  sandy  loam  soil  (Laliberte  et  al., 
1966)  for  which  the  deterministic  and  statistical  properties  are  presented  in  Table  2.4 
and  Table  2.3,  respectively.  From  Gelhar  (1993,  p.  291),  an  averaged  value  of  1  was 
set  for  the  variance  of  InfC*,  afnK,  =  /  =  1.  Uniform  soils  simulations  were  conducted 
by  setting  /  =  0  and  K,*G  =  JC*. 

The  rainfall  series  were  generated  using  WGEN,  a  weather  generator  developed 
by  Richardson  and  Wright  (1984).  This  computer  simulation  code  generates  daily 
rainfall  series  which  reproduces  the  statistics  of  the  rainfall  from  a  specified  region. 
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WGEN  uses  a  first-order  Markov  process  to  represent  the  relative  occurrences  of 
wet  day  (Pw)  and  dry  days  (PD)  and  a  gamma  distribution  to  represent  the  daily 
volume  of  the  rainfall.  Statistical  parameters  of  the  rainfall  series  of  Tampa  and 
Olympia  and  saturated  hydraulic  conductivity  are  summarized  in  Table  2.4.  In  order 
to  compare  the  analytical  development  with  the  numerical  experiments,  a  negative 
exponential  function  was  fit  to  the  sample  covariance  function  produced  by  WGEN 
for  each  weather  scenario.  Note  that  the  annual  average  of  the  actual  monthly  rainfall 
statistics  were  implemented  to  maintain  stationarity  of  the  rainfall  time  series. 

Values  of  soil- water  content,  Darcy's  flux  and  pore- water  velocity  were  recorded 
at  one  node  located  in  the  middle  of  four  control  planes  CP1,  CP2,  CP3,  and  CP4 
and  relevant  statistics  calculated  for  each  of  the  400  days  and  over  all  the  runs.  The 
four  control  planes  are  located  at  .6125,  2.125,  4.125,  and  6.125  m  from  the  land 
surface  for  both  soils. 

2.4.2    Results  and  Discussion 

The  evolution  of  predicted  (analytical)  means  and  standard  deviations  of  the 
soil-water  content,  Darcy's  flux,  and  pore-water  velocity  over  depth  are  compared  to 
their  measured  (Monte-Carlo)  counterparts  for  uniform  soils  in  Figures  2.14-2.16  and 
for  heterogeneous  soils  in  Figures  2.17-2.19. 

Figure  2.14  indicates  that  the  predictions  of  the  stochastic  unsaturated  flow 
model  are  good  in  the  uniform  sandy  loam  soil  subject  to  the  Olympia  weather  pat- 
tern. The  mean  and  the  standard  deviation  of  the  predicted  and  measured  soil-water 
content,  Darcy's  flux  and  pore- water  velocity  are  very  similar  both  at  shallow  depths 
and  at  deep  depths.  In  the  uniform  fine  sandy  soil  subject  to  Olympia  rainfall  pattern 
(Figure  2.15),  the  predictions  remain  reasonably  good.  The  means  of  the  soil-water 
content,  Darcy's  flux  and  pore-water  velocity  are  well  predicted.  The  predictions  of 
the  standard  deviation  are  also  quite  good,  however  the  Darcy's  flux  and  pore-water 
velocity  standard  deviation  do  not  match  quite  as  well  as  for  the  sandy  loam  soil.  The 
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model  performance  deteriorates  for  the  fine  sandy  soil  subject  to  the  Tampa  rainfall 
pattern  (Figure  2.16).  Although  the  means  predictions  of  the  soil-water  content, 
Darcy's  flux  and  pore-water  velocity  remain  of  good  quality,  the  standard  deviations 
tend  to  be  overpredicted  compared  to  the  Monte-Carlo  experiments,  particularly  for 
the  Darcy's  flux  and  pore-water  velocity. 

From  such  results,  it  can  be  concluded  that  the  stochastic  unsaturated  flow 
model  gives  the  best  predictions  in  fine  textured  soils  subject  to  a  rainfall  pattern 
with  small  variance.  In  fact,  the  soil  type  as  well  as  the  rainfall  pattern  explain 
these  observations.  For  a  given  rainfall  pattern,  the  larger  buffering  effect  of  the  finer 
textured  soil  on  the  flow  variability  is  anticipated.  Finer  textured  soils  have  slower 
drainage  which  yields  a  slower  redistribution  process  that  dampens  further  the  pre- 
cipitation events  and  reduces  the  temporal  variability  of  the  soil-water  content.  The 
result  is  that  temporal  fluctuations  around  the  mean  soil- water  content  are  smaller, 
resulting  in  an  inhibition  of  the  nonlinear  behavior  of  the  soil-moisture  diffusivity  and 
hydraulic  conductivity,  and  an  increase  of  the  accuracy  of  the  linearization  of  these 
functions.  For  a  given  soil  type,  the  model  accuracy  also  depends  on  the  rainfall  pat- 
tern. A  rainfall  pattern  characterized  by  a  high  variance  and  long  periods  of  dry  days 
(Tampa)  enhances  the  temporal  variability  of  the  soil- water  content,  particularly  at 
shallow  depths.  Thus,  it  increases  the  nonlinear  behavior  of  the  soil-moisture  diffu- 
sivity and  hydraulic  conductivity.  Thus,  the  best  results  are  obtained  with  the  sandy 
loam  soil  subject  to  Olympia  rainfall  pattern  and  the  worst  results  are  obtained  in 
the  fine  sandy  soil  subject  to  the  Tampa  rainfall  pattern. 

Figures  2.17-2.19  compare  the  measured  statistics  of  the  soil-water  content, 
Darcy's  flux,  and  pore-water  velocity  to  the  predicted  statistics  of  the  stochastic 
unsaturated  flow  model  in  heterogeneous  soils.  These  figures  show  that  the  predicted 
statistics  generally  compare  favorably  to  their  numerical  Monte-Carlo  counterparts. 
The  general  trend  is  faithfully  reproduced  for  the  soil- water  content,  Darcy's  flux  and 
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pore-water  velocity  statistics.  The  quality  of  the  fit  is  better  in  coarse  textured  soils 
than  in  fine  textured  soils,  as  indicated  in  Figures  2.17  and  2.18,  and  improves  as 
the  variance  of  the  rainfall  decreases.  A  fine  textured  soil  tends  to  enhance  lateral 
flow  and  thereby  may  create  two-dimensional  flow  that  is  not  captured  by  the  one 
dimensional  stochastic  unsaturated  flow  model.  As  with  the  case  of  homogeneous 
soil,  the  better  performance  is  observed  for  the  lower  variance  of  the  Olympia  rainfall 
scenario  compared  to  that  of  the  Tampa  rainfall  pattern. 

The  effect  of  the  soil  texture  and  rainfall  pattern  on  the  flow  statistics  explained 
previously  under  uniform  soil  conditions  also  apply  to  the  results  reported  under  het- 
erogeneous soil  conditions.  However,  the  addition  of  the  soil  variability  complicates 
the  flow  pattern  and  creates  a  new  source  of  modeling  error.  In  a  homogeneous  soil 
subject  to  spatially  uniform  but  temporally  random  rainfall,  the  horizontal  pressure 
head  gradient  is  zero  leading  to  no  horizontal  component  of  the  Darcy's  flux.  On 
the  contrary,  in  a  heterogeneous  soil  the  possibility  of  horizontal  gradients  can  not 
be  overlooked.  Soil  heterogeneity  contributes  to  a  non-zero  horizontal  Darcy's  flux 
that  is  smaller  to  that  of  the  vertical  flux,  but  that  can  not  be  neglected.  As  a  con- 
sequence, the  flow  pattern  in  a  heterogeneous  soil  subject  to  spatially  uniform  and 
temporally  random  rainfall  is  more  likely  to  be  two-dimensional.  The  investigation 
made  by  Russo  et  al.  (1994a)  supports  this  flow  analysis.  Running  a  single  numerical 
experiment  of  transport  in  a  heterogeneous  soil,  this  study  found  that  under  non- 
continuous  infiltration  due  to  periodic  spatially  uniform  influx,  the  flow  pattern  was 
clearly  two-dimensional. 

A  graphical  comparison  between  the  Monte-Carlo  simulations  performed  with 
the  homogeneous  soils  (Figures  2.14-2.16)  and  those  performed  with  soil  heterogeneity 
(Figures  2.17-2.19)  illustrate  the  fact  that  i)  in  both  cases  the  measured  flow  statistics 
generally  decrease  over  depth  and  ii)  they  are  slightly  larger  and  more  unstable  in 
the  case  of  soil  and  rainfall  heterogeneity.  The  similarity  of  the  measured  moments 
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between  the  two  cases  underscores  the  dominance  of  the  rainfall  variability  over  that 
of  the  soil  to  the  overall  uncertainty.  Such  results  agree  with  the  study  of  Jury 
and  Gruber  (1989)  that  investigated  the  individual  effect  of  rainfall  variability  and 
soil  heterogeneity  on  pesticide  leaching  potentials  by  coupling  the  climatic  variability 
model  of  Eagleson  (1978)  with  the  soil  variability  transport  model  of  Jury  (1982). 
They  found  that  the  rainfall  plays  a  larger  role  than  soil  variability  in  controlling  the 
movement  of  mobile  and  persistent  pesticides  through  the  vadose  zone. 

Additional  explanations  can  shed  light  on  such  results.  The  developed  stochas- 
tic model  is  based  on  a  one-dimensional  flow  equation  and  requires  linearization  of 
the  soil-moisture  diffusivity  and  soil  hydraulic  conductivity  around  the  mean  soil- 
water  content.  The  Monte-Carlo  simulations  solve  the  two  dimensional  Richards 
equation  which  acknowledges  the  full  nonlinearity  of  the  soil  properties.  Therefore, 
it  is  important  to  note  that  it  is  difficult  to  distinguish  the  discrepancies  between  the 
predicted  and  measured  ensemble  moments  due  to  soil  heterogeneity  from  that  due 
to  the  linearization.  Furthermore,  the  fixed  water  table  at  10  meters  in  the  numerical 
experiment  contrasts  with  the  semi-infinite  domain  described  in  equations  (2.3)-(2.5). 
Another  possible  source  of  error  is  the  fact  that  the  covariance  of  the  discrete  rainfall 
series  generated  by  WGEN  is  approximated  by  a  continuous  distribution.  The  stan- 
dard deviation  predictions  have  been  performed  at  the  second  order.  Previous  works 
have  shown  the  limitations  of  this  low  order  approximation  when  the  input  data  ex- 
hibit coefficient  of  variation  greater  that  the  unity.  Undoubtely,  similar  results  are 
reported  in  the  current  work.  The  evaluation  of  the  statistics  to  second  order  implies 
that  only  the  first  and  second  moments  of  the  rainfall  distribution  are  important  in  de- 
termining the  flow  statistics  in  the  unsaturated  zone.  Finally,  although  the  flow  mass 
balances  of  the  numerical  experiments  never  exceeded  10%,  it  should  be  noted  that 
the  mass  balance  statistics  compare  two  approximate  quantities:  the  time-integrated 
soil-water  content  flux  accross  all  boundaries  and  the  total  mass  change  in  soil-water 
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content  in  the  computational  domain.  When  the  rainfall  pattern  are  highly  variable, 
these  quantities  are  difficult  to  compute,  and  their  difference  may  not  be  a  reliable 
measure  of  numerical  accuracy.  The  measured  means  Darcy's  flux  in  Figures  2.15b 
and  2.16b  illustrate  the  numerical  error  of  the  Monte-Carlo  runs.  The  mean  Darcy's 
flux  from  the  Monte-Carlo  simulation  is  not  exactly  equivalent  to  the  mean  rainfall 
flux  throughout  the  profile,  which  is  an  expected  result  coming  from  the  uniformity 
of  soil  properties. 

2.5  Summary 

The  purpose  of  this  Chapter  was  to  develop  and  test  a  semi-analytical,  one- 
dimensional,  stochastic  model  for  transient  water  flow  in  a  partially  saturated  and 
heterogeneous  soil.  The  model  predicts  the  mean  and  covariance  of  the  soil- water 
content,  Darcy's  flux,  and  pore- water  velocity  as  a  function  of  the  rainfall  and  satu- 
rated hydraulic  conductivity  statistics.  The  one-dimensional  soil- water  content-based 
Richards  equation  for  the  flow  moments  is  solved  in  the  Laplace  domain  after  a  lin- 
earization of  the  relative  hydraulic  conductivity  and  soil-moisture  diffusivity  about  the 
space-invariant  mean  soil- water  content.  Then,  the  analytical  results  are  numerically 
inverted  using  a  modified  Fast  Fourier  Transform  algorithm.  The  model  predictions 
were  compared  to  results  obtained  from  intensive  two-dimensional  Monte-Carlo  sim- 
ulations for  two  different  rainfall  patterns  and  two  different  soil  types. 

The  main  finding  of  this  chaper  are  summarized  below: 

-  At  shallow  depths  in  the  vadose  zone,  both  the  semi-analytical  results  and  the 
Monte-Carlo  simulations  show  that  the  statistics  of  the  soil-water  content,  Darcy's 
flux  and  pore-water  velocity  are  dominated  by  the  effect  of  rainfall  variability  rather 
than  by  the  effect  of  the  soil  heterogeneity.  At  deeper  depths  where  the  temporal 
variability  of  the  rainfall  is  attenuated,  the  soil  heterogeneity  becomes  more  dominant. 
-  Compared  to  the  Monte-Carlo  simulations,  the  performance  of  the  stochastic  flow 


53 


model  can  be  evaluated  by  comparing  two  extreme  cases.  The  statistics  of  the  soil- 
water  content,  Darcy's  flux  and  pore- water  velocity  are  much  better  predicted  in  a  fine 
textured  soil  subject  to  a  small  rainfall  variance  than  in  coarse  textured  soil  subject 
to  a  high  rainfall  variance.  The  fundamental  model  assumption  of  a  small  variability 
of  the  soil-water  content  around  a  constant  mean  value  explains  such  results.  Any 
driving  force  that  increases  significantly  the  variability  of  the  soil-water  content  and 
thus  the  non-linear  of  the  flow  system  will  tend  to  violate  this  assumption.  The 
rainfall  variability  is  the  most  critical  driving  force  that  affects  the  temporal  evolution 
of  the  soil-water  content  at  shallow  depths  and  therefore  the  nonlinearity  of  the  flow 
equation.  To  a  lesser  extent,  the  storage  capacity  of  soil  is  a  secondary  driving 
force  in  the  sense  that  it  will  attenuate  the  effect  of  the  rainfall.  Indeed,  the  better 
predictions  of  flow  statistics  are  observed  with  the  homogeneous  sandy  loam  soil 
than  with  the  homogeneous  fine  sandy  soil  because  the  higher  storage  capacity  and 
slower  redistribution  tend  to  reduce  the  temporal  variability  of  the  rainfall.  Note  that 
because  the  stochastic  flow  model  is  one  dimensional,  its  performances  are  better  in 
coarse  textured  soils  than  in  fine  textured  soils  for  the  heterogeneous  cases. 

-  The  inversion  of  the  double  Laplace  transform  using  a  modified  Fast  Fourier 
Transform  algorithm  is  a  quick  and  accurate  method  to  estimate  the  moments  of  soil- 
water  content,  Darcy's  flux  and  pore- water  velocity  that  is  much  more  computationaly 
efficient  than  the  Monte-Carlo  simulations.  It  can  be  roughly  estimated  that  comput- 
ing the  flow  statistics  for  each  Monte-Carlo  run  generally  needs  iV2  order  operations 
per  time  step  for  N  discrete  nodes  for  a  two-dimensional  problem.  The  numerical  Fast 
Fourier  Transform  involves  only  on  the  order  of  4  *  Log2N  operations  for  N  nodes. 
The  CPU  time  used  for  these  two-dimensional  randomly  transient  boundary  prob- 
lems on  a  Sun  Ultra  workstation  is  about  3  minutes  for  the  Fast  Fourier  Transform 
approximation  versus  twenty-five  days  using  a  set  of  150  Monte-Carlo  runs. 
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Some  changes  can  be  made  to  improve  the  performance  of  this  model  as  well  as 
to  extend  its  capability.  The  systematic  approach  based  on  the  converging  series  offers 
the  advantage  of  tracking  the  order  of  estimation.  Computing  ensemble  statistics  of 
known  order  is  a  consequence  of  this  method.  Thus  improvement  of  estimates  are 
possible  by  considering  higher  orders.  However,  such  endeavor  would  require  solving 
more  complex  equations  that  may  rely  exclusively  on  numerical  methods.  Another 
possible  improvement  would  be  to  extend  this  model  to  two  and  three  dimensions,  to 
integrate  a  spatio-temporally  variable  rainfall  series  and  to  incorporate  water  uptake 
in  the  governing  flow  equation.  This  latter  improvement  is  undertaken  in  the  next 
Chapter,  and  the  former  task  is  left  for  future  work. 

Several  applications  exist  for  this  developed  model.  It  is  becoming  increasingly 
recognized  that  the  soil-water  content  plays  a  determinant  role  in  predicting  global 
climatic  patterns.  Soil  water  content  affects  the  vertical  fluxes  of  energy  and  moisture 
to  the  atmosphere  and  deeper  subsurface  and  also  the  horizontal  fluxes  of  moisture, 
namely  runoff  {Budyko,  1956,  Delworth  and  Manabe,  1988).  Knowledge  of  the  soil- 
water  content  statistics  should  improve  the  parameterization  of  land-surface  modeling. 
Finally,  the  statistics  of  the  pore- water  velocity  can  be  used  to  predict  solute  transport 
in  soil.  This  topic  is  addressed  in  Chapter  4. 
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Table  2.1:  Deterministic  hydraulic  properties  of  the  fine  sandy  soil. 


Property 

Value 

es 

0.35 

eT 

0.063 

K 

-.82  m 

X 

3.7 

Table  2.2:  Statistical  properties  of  the  rainfall  series  and  of  the  Log-saturated  hy- 
draulic conductivity. 


Statistics- Variable 

Rainfall 

Ln  k; 

Mean 
Variance 

Correlation  length 

0.0038  m/d 
9.81  x  10_5m2/d2 
0.6  d 

1.61 
1 

1  m 

Table  2.3:  Deterministic  hydraulic  properties  of  the  Brooks-Corey  model  for  the  fine 
sandy  soil  and  sandy  loam  soils. 


Property 

Fine  sand 

Sandy  loam 

o. 

.35 

.496 

8r 

.063 

.11 

X 

3.7 

1.6 

h* 

-.82  m 

-.85  m 
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Table  2.4:  Ensemble  statistics  of  the  two  rainfall  scenarios  combined  with  the  two 
soils,  SI:  Tampa  rainfall  x  fine  sandy  soil,  S2:  Olympia  rainfall  x  fine  sandy  soil,  S3: 
Olympia  rainfall  x  sandy  loam  soil.  Pi[W/D]  and  P^W/W]  stand  for  the  probabilities 
of  a  wet  day  on  day  i  given  a  dry  day  on  day  i-1,  and  the  probability  of  a  wet  day 
on  day  i  given  a  wet  day  on  day  i-1,  respectively.  o^.,  C%5,  and  Xs  stand  for 
the  mean,  the  variance,  the  coefficient  of  variation  and  the  temporal  correlation  scale 
of  the  rainfall  series,  respectively,  and  JC*G,  a?njc<,  and  A„  represent  the  mean,  the 
variance,  and  the  spatial  correlation  scale  of  the  log-saturated  hydraulic  conductivity, 
respectively. 


SI 

S2 

S3 

Pi[W/W] 

.453 

.816 

.816 

Pi[W/D] 

.240 

.452 

.452 

Pw 

.30 

.71 

.71 

fag)  (m/d) 

.00384 

.00293 

.00293 

a\.  (xlO-5m'V^) 

9.81 

1.76 

1.76 

2.58 

1.43 

1.43 

X5  (d) 

.6 

.6 

.6 

K*G  (m/d) 

5. 

5. 

.7 

1 

1 

1 

Xn  (m) 

1 

1 

1 

#00 

.1008 

.0980 

.216 
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Normalized  soil-water  content  autocovariance 


Figure  2.1:  The  normalized  soil- water  content  autocovariance,  Pm^z  f™^2  '^m't\  at 

°  '  <Tg(Z  =  Z  —  0)  ' 

z=z'=0.5  meter. 
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0.0  0.2  0.4  0.6  0.8  1.0 
Normalized  soil-water  content  autocovariance 


Figure  2.2:  The  normalized  soil-water  content  autocovariance  at  z=z'=0.5  meter, 
pBe(z-.5mj£-^5m,t ) ^  ^  ^e  absence  0f  rainfall  variability. 
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i  i     1     i     1     i     1     I     '  i 
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Normalized  soil-water  content  autocovariance 


Figure  2.3:  The  normalized  soil-water  content  autocovariance,  pw(z-4™^f -*">■<*')  a^ 

ag\z— 2  —  "I 

z=z'=4  meters. 
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0.0       0.2       0.4       0.6       0.8  1.0 
Normalized  Darcy's  flux  autocovariance 


Figure  2.4:  The  normalized  Darcy's  flux  autocovariance,  Pn{z~^=$'^jm,e) ,  at 
z=z'=0.5  meter. 
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0.0       0.2        0.4        0.6       0.8  1.0 
Normalized  pore-water  velocity  autocovariance 


Figure  2.5:  The  normalized  pore-water  velocity  autocovariance,  Pvv^z~i™*%'~^m,1t\ 
at  z=z'=0.5  meter. 
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1.6 

2.1  z',  meters 
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Normalized  soil-water  content  autoco variance 


Figure  2.6:  The  normalized  soil- water  content  autocovariance,  pgg(*'*-224*V'-22<*)  at 
t=t'=22  days. 
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Normalized  soil-water  content  autocovariance 


Figure  2.7:  The  normalized  soil-water  content  autocovariance,  pM^t-^d-,z\t'-22d)  at 
t=t'=22  days  in  the  absence  of  rainfall  variability. 
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Figure  2.8:  The  normalized  soil- water  content  autocovariance,  p»»iz<t-^od'f't'-60d)  at 

<Xj(z=z'=0)  ' 

t=t'=60  days. 
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Normalized  Darcy's  flux  autocovariance 


Figure  2.9:  The  normalized  Darcys's  flux  autocovariance,  pw(z+-^*'^-22d)  ;  at 
t=t'=22  days. 
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Normalized  pore-water  velocity  autocovariance 


Figure  2.10:  The  normalized  pore- water  velocity  autocovariance,  F""^2(^f='o)  ^  • 
at  t=t'=22  days. 
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Figure  2.11:  a)  Effect  of  the  variance  and  of  the  correlation  scale  of  the  saturated 
hydraulic  conductivity  on  the  temporal  evolution  of  the  soil-water  content  autocor- 
relation at  z=z'=0.4  meter,  b)  Effect  of  the  variance  and  of  the  correlation  scale  of 
the  saturated  hydraulic  conductivity  on  the  spatial  evolution  of  the  soil-water  content 
autocorrelation  at  t=t'=  60  days. 
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Figure  2.12:  a)  Effect  of  the  variance  and  of  the  correlation  scale  of  the  saturated 
hydraulic  conductivity  on  the  temporal  evolution  of  the  pore-water  velocity  autocor- 
relation at  z=z'=0.4  meter,  b)  Effect  of  the  variance  and  of  the  correlation  scale  of  the 
saturated  hydraulic  conductivity  on  the  spatial  evolution  of  the  pore-water  velocity 
autocorrelation  at  t=t'=  60  days. 
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Figure  2.13:  Physical  description  of  the  infiltration  experiment 
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Figure  2.14:  Predicted  versus  measured  mean  and  standard  deviation  of  the  a)  soil- 
water  content,  b)  Darcy's  flux  and  c)  pore- water  velocity  in  the  middle  of  each  control 
plane  for  the  uniform  sandy  loam  soil  combined  with  the  Olympia  rainfall  scenario. 
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Figure  2.15:  Predicted  versus  measured  mean  and  standard  deviation  of  the  a)  soil- 
water  content,  b)  Darcy's  flux  and  c)  pore-water  velocity  in  the  middle  of  each  control 
plane  for  the  uniform  fine  sandy  soil  combined  with  the  Olympia  rainfall  scenario. 
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Figure  2.16:  Predicted  versus  measured  mean  and  standard  deviation  of  the  a)  soil- 
water  content,  b)  Darcy's  flux  and  c)  pore-water  velocity  in  the  middle  of  each  control 
plane  for  the  uniform  fine  sandy  soil  combined  with  the  Tampa  rainfall  scenario. 
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Figure  2.17:  Predicted  versus  measured  mean  and  standard  deviation  of  the  a)  soil- 
water  content,  b)  Darcy's  flux  and  c)  pore-water  velocity  in  the  middle  of  each  control 
plane  for  the  heterogeneous  sandy  loam  soil  combined  with  the  Olympia  rainfall 
scenario. 
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Figure  2.18:  Predicted  versus  measured  mean  and  standard  deviation  of  the  a)  soil- 
water  content,  b)  Darcy's  flux  and  c)  pore- water  velocity  in  the  middle  of  each  control 
plane  for  the  heterogeneous  fine  sandy  soil  combined  with  the  Olympia  rainfall  sce- 
nario. 
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Figure  2.19:  Predicted  versus  measured  mean  and  standard  deviation  of  the  a)  soil- 
water  content,  b)  Darcy's  flux  and  c)  pore-water  velocity  in  the  middle  of  each  control 
plane  for  the  heterogeneous  fine  sandy  soil  combined  with  the  Tampa  rainfall  scenario. 


CHAPTER  3 

UNSATURATED  TRANSIENT  STOCHASTIC  FLOW  MODEL  WITH  WATER  UPTAKE 

3.1  Introduction 

Realistic  description  of  unsaturated  flow  in  soil  must  consider  the  important 
role  played  by  water  uptake.  Through  its  two  forms,  evaporation  and  transpiration, 
water  uptake  modifies  the  flow  pattern  over  time  and  space  during  redistribution  and 
depletes  the  soil-water  content.  For  example,  evapotranspiration  averages  out  to  be 
87%  of  the  annual  precipitation  in  Australia  (Australian  Water  Resources  Council, 
1976  )  and  72%  in  Florida  (Jones  et  al,  1984)-  As  a  result,  water  uptake  diminishes 
the  level  of  recharge  and  also  delays  or  reduces  the  transport  of  contaminants  to 
groundwater. 

Few  studies  accounting  for  the  effect  of  water  uptake  on  the  flow  statistics  in 
the  unsaturated  zone  have  been  reported.  Several  factors  contribute  to  this  observa- 
tion. Water  uptake  has  been  assumed  not  to  vary  over  space  as  much  as  saturated 
hydraulic  conductivity,  and  the  random  temporal  changes  of  water  uptake  have  been 
assumed  to  be  less  important  than  those  of  precipitation.  Thus  in  most  reported 
studies  on  spatio-temporal  variability  of  water  flow  in  the  vadose  zone,  water  uptake 
has  been  neglected.  However,  it  is  well  established  that  evaporation  and  transpiration 
are  function  of  soil-water  content,  which  depends  upon  the  hydraulic  conductivity. 
Therefore,  spatial  variability  of  water  uptake  is  indirectly  related  to  saturated  hy- 
draulic conductivity  variability  through  soil-water  content.  Some  previous  works  have 
accounted  for  the  water  uptake  in  the  prediction  of  water  flow.  Warrick  (1974)  devel- 
oped deterministic  analytical  solutions  of  the  steady-state  one-dimensional  linearized 
unsaturated  flow  equation  with  a  spatially  variable  water  uptake  function  expressed  as 

76 


77 


a  negative  exponential  function  of  depth.  Rubin  and  Or  (1993)  developed  and  tested 
a  one-dimensional,  steady-state,  stochastic-analytical  flow  model  that  also  included 
an  exponentially  decaying  water  uptake  function  over  depth.  Their  analysis  consisted 
of  examining  the  effect  of  a  spatially  random  water  uptake  depth  on  predicting  the 
variances  of  pressure  head  and  saturation  level  in  a  heterogeneous  soil  overlying  a 
water  table.  They  found  that  water  uptake  causes  a  sharp  decrease  in  the  variances 
of  the  pressure  and  saturation  level  from  their  stationary  value  near  the  surface  to 
zero  at  the  water  table. 

Katul  et  al.  (1997)  developed  a  transient,  first-order  water  budget  equation 
for  a  volume- averaged  soil-water  content  that  included  a  random  root  water  uptake 
expressed  as  a  function  of  the  soil-water  content  and  a  deterministic  potential  transpi- 
ration rate.  They  solved  the  equation  for  the  time-space  soil-water  content  variance 
and  tested  it  against  some  field  measurements  and  numerical  simulations  where  the 
soil  properties  were  random.  Their  analysis  showed  that  two  dynamic  processes  are 
responsible  for  the  time  variation  of  the  soil-water  content  spatial  variance  i)  the  in- 
teraction between  the  soil-water  content  and  the  vertical  flux  gradient,  and  ii)  the 
interaction  between  the  soil-water  content  and  the  root  uptake.  They  also  found  that 
the  root  uptake  component  of  the  variance  is  comparable  to  the  contribution  from 
soil  properties  and  that  soil  water  redistribution  and  is  central  to  the  dissipation  of 
the  soil- water  content  variance. 

More  recently,  Russo  et  al.  (1998)  performed  a  single  numerical  three-dimensional 
simulation  of  unsaturated  flow  and  transport  under  a  transient  deterministic  weather 
boundary  condition  and  a  water  uptake  flux  throughout  the  soil  profile.  The  root 
water  uptake  was  a  function  of  the  local  soil-water  content,  a  deterministic,  tempo- 
rally varying  transpiration  rate,  and  a  deterministic  root  effectiveness  function  which 
changed  over  depth  but  was  time-invariant.  They  also  assumed  that  the  saturated 
hydraulic  conductivity  was  spatially  variable.  Among  their  findings,  their  analysis 
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showed  that  the  water  uptake  increases  the  spatial  variances  of  the  soil-water  content 
and  longitudinal  velocity,  and  causes  steeper  head  gradients. 

These  studies  discussed  above  are  restricted  by  severe  simplifying  assumptions 
regarding  the  unsaturated  flow  regime  (very  often  assumed  to  be  steady  state  or  de- 
terministic), by  the  limited  number  of  numerical  simulations  and  by  their  limited 
scope  regarding  the  flow  parameters  under  investigation.  In  this  Chapter,  considera- 
tion of  temporally  random  rainfall  and  water  uptake  in  association  with  the  spatial 
variability  of  the  saturated  hydraulic  conductivity  will  be  accounted  in  predicting  the 
statistics  of  soil-water  content,  Darcy's  flux  and  pore-water  velocity.  The  objectives 
of  this  Chapter  are  to  i)  derive  a  one-dimensional  stochastic  flow  model  in  heteroge- 
neous soil  subject  to  transient  flow  regime  due  to  rainfall  and  water  uptake,  and  ii) 
to  examine  the  effect  of  a  spatially  variable  and  temporally  random  water  uptake  on 
the  statistics  of  the  soil-water  content,  Darcy's  flux  and  pore-water  velocity. 

3.2    Analytical  Development 

In  soils,  rainfall  and  irrigation  patterns  as  well  as  water  uptake  are  the  driving 
forces  for  water  flow.  Gravitationally  dominated  flow  is  often  assumed  to  take  place 
in  unsaturated  soils.  This  simplification  of  a  unidimensional  flow  has  been  frequently 
reported  in  coarse  textured  soils  and  with  a  high  infiltration  level  compared  to  the 
saturated  hydraulic  conductivity  (Butters  et  ai,  1989;  Russo,  1991 ;  Katul  et  ai, 
1991).  Decreasing  quantity  of  extracted  water  over  depth  is  often  observed  in  natural 
plant  systems.  This  is  why  water  uptake  has  been  previously  modeled  under  a  steady 
state  flow  regime  using  a  negative  exponential  weighting  function  over  depth  (Raats, 
1974;  Warrick,  1974;  Rubin  and  Or,  1993).  Taking  advantage  of  these  modeling 
simplifications,  a  one-dimensional  equation  for  transient  unsaturated  water  flow  with 
a  water  sink  term  is  assumed.  The  saturated  hydraulic  conductivity  and  water  uptake 
are  assumed  to  be  randomly  variable  in  the  governing  equation,  and  rainfall  is  assumed 
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to  be  randomly  variable  in  the  surface  boundary  condition.  Thus,  in  a  rigid  soil 
matrix,  the  one-dimensional  dimensionless  water  flow  equation  in  variably  saturated 
and  heterogeneous  conditions  subject  to  transient  flow  regime  due  to  rainfall  and 
water  uptake  is  given  by: 


where  JCs(z)  T>(9)  is  the  soil-moisture  diffusivity,  JC(6)  is  the  relative  hydraulic 
conductivity  which  is  a  function  of  the  local  volumetric  soil- water  content  9,  ICs(z) 
is  a  spatially  random  saturated  hydraulic  conductivity,  £  is  the  mean  water  uptake 
rate  over  the  rooting  depth,  zc  is  the  characteristic  rooting  depth,  £  is  a  small  scaling 
parameter,  and  x(t)  is  a  random  function  of  time  with  a  zero  mean  and  unit  standard 
deviation.  Note  that  the  Darcy's  flux  is  oriented  positively  downward  along  the 
vertical  direction  z  and  is  defined  as: 
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q  =  -Ka{z)V{6)—  +  Ka{z)K{9)  (3.2) 

In  order  to  account  for  soil  variability,  a  lognormal  distribution  was  assumed  for 
the  saturated  hydraulic  conductivity,  JCs*(z)  =  K,*Ge^^  where  K,q*  is  the  geometric 
mean  [LT_1],  /  is  a  constant  parameter  equivalent  to  the  standard  deviation  of  In 
K*  and  rj  is  a  normally  distributed  random  space  function  of  zero  mean  and  unit 
variance.  It  is  further  assumed  that  all  these  quantities  are  relevant  at  the  scale 
of  the  representative  elementary  volume  and  that  the  air  phase  is  immobile.  The 
previous  equations  (3.1)  and  (3.2)  are  made  dimensionless  by  using  the  following 
scaling  factors: 
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where  /i£  denotes  the  local  bubbling  pressure  [L],  (e*)  and  ae*  are  the  dimen- 
sional mean  and  standard  deviation  of  the  water  uptake  [LT'1],  respectively.  The 
asterisks  indicate  the  dimensional  quantities,  and  the  ()  stands  for  the  expected  value 
operator.  Note  also  that  the  current  water  uptake  model  differs  from  that  used  by 
Raats  (1974),  Warrick  (1974)  and  Rubin  and  Or  (1993)  in  the  sense  that  it  is  the 
product  of  a  negatively  exponential  depth  weighting  function  by  a  temporally  random 
flux  decomposed  in  terms  of  a  stationary  mean  and  a  temporally  random  perturba- 
tion. At  the  bottom  boundary  condition,  the  Darcy's  flux  is  deterministic  and  its 
value  is  set  to  be  equal  to  the  average  net  applied  water. 

Ignoring  surface  runoff  and  assuming  a  deep  water  table,  the  dimensionless 
initial  and  boundary  conditions  for  a  semi-infinite  domain  are: 


-K,{z)V{d)^-z+K,{z)1C{0) 


=  ll  +  e5{t)  (3.3) 

z=0,t 


BO 


=  n-S  (3.4) 

z=oo,t 

0lt=0   =  0i  (3-5) 


where  TZ  is  the  dimensionless  mean  rainfall  defined  as  e  is  defined  as 

G 

with  aq'Q  representing  the  standard  deviation  of  the  rainfall  series,  and  S(t)  is  a 
random  function  of  time  with  a  zero  mean  and  unit  standard  deviation. 

3.2.1    Power  Scries  Expansion 

The  difficulty  of  solving  the  nonhomogeneous  partial  differential  equation  gov- 
erning flow  in  unsaturated  porous  media  rests  in  the  nonlinear  coefficients.  We  assume 


81 


that  the  solution  to  the  nonlinear  partial  differenctial  equation  can  be  approximated 
by  a  power  series  solution  of  the  form  (Met,  1997): 


oo      oo  oo 


6(z,  W,C)  =  EEE  0ljk(z,  t)  eTC*  (3-6) 

i=0  j=0  k=0 

where  the  equation  is  non-dimensionalized  so  that  e,  /,  and  (  are  small  di- 
mensionless  parameters  and  thus  convergence  of  the  series  is  insured.  Following  the 
procedure  outlined  in  the  previous  Chapter,  the  non-linear  functions  V{9)  and  1C(9) 
are  expanded  about  #0oo  in  a  Taylor's  series  as  follows: 


g(0)  =  g(e0oo)  +  % 


(0-0ooo)+^ 

"ooo 


(9  -  Oqoo  ) 2 


where  each  of  the  partial  derivatives  is  to  be  evaluated  at  0Ooo-  Note  again 
that  each  term  (6  —  #ooo)n  can  be  simplified  as: 

(00      oo      oo  \  71         /  oo      oo  oo 

E  E  E  M*.  *)  w  -  ^ooo  =  E  E  E M*. *) 
i=0  j=0  k=0  /  \i=l  j=l  k-1 

Substituting  this  power  series  expansion  in  equations  (3.1)  and  (3.3)-(3.5),  the 
left  hand  side  in  equation  (3.1)  becomes  to  the  second  order: 


p.,  9ooo  +  e  #ioo  +  /  #oio  +  C^ooi  +  e/  0no  +     #101  +  \ 

30  _     1  K  0QU  +  <?  0200  +  P  0020  +  C20QO2  +  •  ■  ■  ' 

dt  ~  at 

(3.7) 


and  the  right  hand  side  is 
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d_ 

dz 


(1  +  fv+ 


V 


/ 


/  V(6000)  +  §\6oqo  (e9m  +  f60W  +  C0OO1+  \ 

e/0no  +  <  0ioi  +  /C  0on+ 
e2  02oo  +  P  0O2o  +  C20oo2  +  •  •  • ) 

(e  ^100  +  /  ^010  +  C^001  + 

e/ 0no +  <  0ioi +/C  0oii  + 

0200  +  /2  0020  +  C20OO2  +  ■  ■  •  )2  +  •  ■  • 

^(0ooo  +  e  0ioo  +  /  0oio  +  C0ooi  + 

e/0no +  <  0ioi  +/C  0oii  + 
e2  0200  +  /2  0020  +  C20oo2  +  ■  •  •  ) 
/  /C(0Ooo)  +  % \9om  (e0ioo  +  /0oio  +  C0ooi+  \ 
ef  0iio  +  <  0ioi  +  /C  0on  + 

0200  +  /2  0020  +  C  0002  +  ■  •  ■  ) 
+  2  fiF  0100  +  / 0010  +  C0OO1+ 

#000 

ef  0iio  +  <  0ioi  +  /C  0oii  + 
e2  0200  +  P  0020  +  C20oo2  +  ---)2  +  -- 

z 


V 


/ 


e  zc 


(£  +  Cx(t) ) 


(3.8) 


where  the  exponential  form  ein  has  been  approximated  by  second  order  Tay- 
lor's series  expansion.  Substituting  the  power  series  expansion  of  6  in  equations 
(3.3)-(3.5),  the  boundary  conditions  become 
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£>(0ooo)  +  wLooo  ^  6,100  +  -^010+  \ 
C0ooi  +  e/0no  +  <0ioi+/C0on  + 
e2  0200  +  f2  0020  +  C20oo2  +  ■  ■  ■  ) 


i        i/lio  i       "101  t  j  s 
*200  +  P  0020  +  C2^002  +  ■  ■  •  ) 

'  2  "Sil^  (e^100  +  /0O1O  + 

#000 

C0ooi+e/0no  +  <0ioi  +/C0on  + 
V  e202oo  +  /20O2o  +  C20oo2  +  ...)2  +  ---  / 
^■(^ooo  +  e^ioo  +  /0O1O  + 
C0ooi+e/0no  +  <0ioi  +  /C0oii  + 
0200  +  P  0020  +  C2^002  +  ...) 

£(0OOo)  +  wLqoo  (6^100  + /^010  + 

C0ooi  +  e/  0no  +  <  0ioi  +  /<  0oii  + 

C2  0200  +  /2  0020  +  C2$002  +  •  ■  ■  ) 
^klxF  (^  #100  +  /  0010  + 

#000 

C#ooi  +  ef  9U0  +  e(  0ioi  +  /C  0on  + 
V  e202oo  +  /20O2o  +  C20oo2  +  ...)2  +  ---  / 


(    P(0ooo)+f  L000  M100  +  / 0oio+  \ 

C^OOl  +6/0110  +  <  0101  +/C  0011  + 
f2  0200  +  /2  0020  +  C20OO2  +  ■  •  •  ) 

^2  ae2 


(e0ioo  +  /0oio+ 

fooo 

C0ooi  +  ef  0iio  +  eC  0ioi  +  K  0on  + 
V  e2  02OO  +  /20o2o  +  C20oo2  +  ..-)2  +  ---  / 

^(0000  +  e0ioo  +  /0O1O  + 

C0ooi  +  ef  0iio  +  <  0ioi  +/C  0oii  + 
e2  0200  +  /2  0020  +  C20oo2  +  •  •  •  ) 
<    A:(0OOo)  +  H|eooo  (e 0ioo  +  / 0oio+  \ 
C0ooi  +  ef  0iio  +  eC  0ioi  +  fC  0oii+ 
e2  0200  +  /2  0020  +  C20oo2  +  •  •  ■  ) 
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^000 
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\  e2  0200  +  /2  0020  +  C20oo2  +  •  •  ■  )2  +  •  •  •  / 


1l  +  eS(z,t) 
at  z  =  0 


(3.9) 


n-s 

at  z  =  00 


(3.10) 


while  the  initial  condition  is 
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#ooo  +  e  #100  +  /  #010  +  C#ooi  +  e/  #110  +     #101  +  fC  #on  + 

#200  +  f2  #020  +  C2#002  +  •  •  • 


) 


z,t=0 


=  0i 


(3-11) 


The  general  solution  of  8(z,  t)  requires  solving  the  governing  flow  equation  and 
its  initial  and  boundary  conditions  for  each  #ijfc,  up  to  the  desired  order. 

3.2.2    Governing  Equations  for  the  Flow  Ensemble  Moments 

The  saturated  hydraulic  conductivity,  the  rainfall  and  the  water  uptake  are 
random  driving  forces  of  the  flow  system.  As  a  result,  the  soil-water  content  becomes 
a  spatio-temporally  random  process.  In  the  absence  of  detailed  measurements  of  these 
forcing  terms,  it  is  impossible  to  characterize  9(z,t).  However,  the  central  tendency 
of  its  behavior  can  be  measured  with  the  ensemble  mean  and  the  expected  variability 
around  that  central  tendency  estimated  using  the  ensemble  variance. 

Soil-water  content  statistics 

Taking  the  expected  value  of  equation  (3.6)  leads  to  the  mean  soil-water  con- 
tent: 


The  zero  order  mean  soil-water  content  is  determined  by  solving  the  follow- 
ing deterministic  equation  that  has  been  established  by  collecting  terms  of  order 


(6(z,  t))  =  (#ooo)  +  e  (#100)  +  /  (#010)  +  C(#ooi)  +  tf  (0UO)  +  <  (#101)  + 


K  (#on)  +  e2  (#200)  +  /2  (#020)  +  C2  (#002)  +  •  •  •  (3.12) 


O(e°/0C°)  in  equations  (3.7)-(3.11): 
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<9#ooo      d  .   .     .  d#ooo 

{~dT ~  d~z[v{6ooo)~dr 

(-v(e00o)^  +  JC(e000) 

(^oooLt=0)  =  @i 


ic(e000)])  =  —5 


2=0,t 


)  =  n-e 


z=+oo,t 


(3.13) 
(3.14) 

(3.15) 
(3.16) 


The  boundary  conditions  considered  here  are  deterministic  and  constant  in 
time  and  create  a  nonlinear  decrease  of  0Ooo  over  depth.  Therefore,  (#ooo)  =  #ooo  and 
away  from  the  initial  condition,  the  previous  equation  simplifies  to  a  nonhomogeneous 
differential  equation  with  nonconstant  coefficients: 


d  \T>(ft    \  dflooo 


—  /C(#ooo)]  =   £ 


— ^(flooo)  d^°Q  +  ^-(^ooo) 
-^(flooo)         +  £(#ooo) 


=  11 


z=0,t 


=  n-s 


(3.17) 
(3.18) 

(3.19) 


2=+00,t 

Once  relative  hydraulic  conductivity  and  soil-moisture  diffusivity  functions 
have  been  specified,  the  zero  order  mean  soil-water  content  can  be  solved  numeri- 
cally. However,  here  the  solution  of  this  equation  is  determined  using  a  simplified 
approach.  The  zero  order  mean  soil- water  content  has  been  estimated  using  the 
equivalent  steady-state  mean  soil  water  (Salvucci  and  Entekhabi,  1994)  defined  as 
#ooo-  This  quantity  corresponds  to  the  mean  soil- water  content  deep  in  the  profile 
where  the  effects  of  rapidly  varying  boundary  conditions  are  damped  out  by  the  over- 
lying soil.  This  mean  soil- water  content  is  therefore  calculated  using  a  quasi-steady 
flow  regime,  assuming  that  the  vertical  gradient  of  the  zero  order  soil-water  content 
is  zero.  Thus,  the  equivalent  steady-state  mean  soil-water  content  is  the  water  con- 
tent that  will  transmit  the  mean  recharge  to  groundwater  (mean  rainfall-mean  water 
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uptake).  This  results  in  a  zero-order  estimate  of  the  mean  soil- water  content  that  is 
space-invariant.  Away  from  the  land  surface,  0OOo  is  calculated  by  solving  equation 
(3.15)  and  noting  that  %*  =  0: 


K{dm)-K  +  S  =  0  (3.20) 

A  consequence  of  assuming  an  equivalent  steady-state  soil- water  content  that  is 
space  invariant  is  that  the  governing  flow  equations  for  higher  order  soil-water  content 
perturbations  become  nonhomogeneous  partial  differential  equations  with  constant 
rather  than  spatially  varying  coefficients.  Therefore,  the  governing  equations  of  01Oo, 
#oio,  #iio,  #2oo,  and  802o  and  their  expected  values  remain  identical  to  these  derived 
in  the  previous  Chapter  for  0iO,  0Oi,  #n,  #20,  and  0O2  and  are  not  repeated  here.  Only 
the  governing  equations  of  the  new  variables  pertaining  to  the  random  water  uptake 
term  0OOi,  #101,  #011,  and  0OO2  are  developed.  Collecting  terms  of  order  0(e°/°C)  in 
equations  (3.7)-(3.11)  leads  to  the  governing  equation  of  900\: 


d6{ 


Tt(  a    ^d20ooi  dK 

dT~v{9ooo)~dz^+  30 
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001 


e  *c 


(-2^000)-^+  a, 

(-V(6000)—+  — 
0ooiL)t=o  =  0 


0ooo 


Sooo 
#001  ) 

0001  ) 


dz 


x(t) 


(3.21) 


=  0 


0000 


z=0,t 


z=+oo,t 


This  governing  equation  resembles  an  advection-dispersion  equation  with  con- 
stant coefficients  and  an  interior,  near  surface  random  forcing  term.  Taking  the 
expected  value  of  equation  (3.21)  leads  to: 
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(^001  )) 


dz 


(x(t)) 


=  0 


z=0,t 


0OOO 


=  0 


2  =+00,t 


Because  the  boundary  and  initial  conditions  are  equivalent  to  zero  and  the 
expected  value  of  the  forcing  term  is  zero,  i.e.  (x)  =  0,  the  expected  value  of  0Ooi 
is  also  zero.  Taking  the  expected  value  after  collecting  terms  of  order  O(e0fC)  in 
equations  (3.7)-(3.11)  leads  to  the  governing  equation  of  (#on): 
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This  governing  equation  resembles  to  a  nonhomogeneous  advection-dispersion 
equation  with  constant  coefficients  and  with  spatio-temporally  varying  boundary  con- 
ditions. Note  that  if  the  cross-covariance  between  saturated  hydraulic  conductivity 
and  water  uptake  is  zero,  then  (#on)  =  0.  Taking  the  expected  value  after  collecting 
terms  of  order  O(ef0C)  in  equations  (3.7)-(3.11)  leads  to  the  governing  equation  of 
(tfioi): 
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Taking  the  expected  value  after  collecting  terms  of  order  C>(e0/°C2)  in  equa- 
tions (3.7)-(3.11)  leads  to  the  governing  equation  of  (#002): 
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where  (flooi^)  simplifies  to  l^f1.  Once  (01Oi),  (0Oii>,  and  (#002)  have  been 
evaluated,  then  higher  order  corrections  to  mean  soil-water  content  can  be  calculated. 
The  autocovariance  of  9,  Pee,  to  the  second  order  is  established  as  follows: 


89 


=  < 


Pee(z,t;z',t')  =((9(z,t)  -  (*(*,*)»  W ,H)  -  (9(z',t')))) 

(       eQOo(z,t)  +  60wo(z,t)  + f6010{z,t)  +  (90oi{z,t)+  \ 
ef  6>no  (2,  t)  +  e(  6Wi  (z,  t)  +  f(  0On  (z,  t)+ 
e2  9200{z,  t)  +  f2  6020(z,  t)  +  (%02(z,  t)  +  . . .  - 
((0ooo(*>  *)>  +  e  (dioo(z,  t))  +  f  (6010(z,  t))  +  C(eQOi(z,  t))+ 
ef  (9lw(z,  t))  +  eC  (eWi(z,  t))  +  K  (6on(z,  t))+ 
V      e2  (02oo(z,  t))  +  f2  (6020(z,  t))  +  C2  <0oo2(2,  *)>  +  ...) 

/  0ooo(2,,O+^ioo(2'^,)  +  /0oio(2',O  +  C0ooi(^,n+  \ 
ef  0no  (z',  t')  +  eC  0ioi  (z',  ?)  +  K  0on  (z',  t')+ 
e2  02oo(^',  f)  +  P  002o(z',  if)  +  C20oo2(^',  0  +  •  •  ■  - 
((e0o0(z',t'))+e(9W0(z',t'))  +  f(6olo(z',t'))+  ) 

C(0ooi(2',  t'))  +  ef  (eno(z',  If))  +  <  (0ioi {z\  ?))  + 
K  (00n(z',  If))  +  e2  (<«*',  f))  +  f2  (0020(2',  t'))+ 
V  C2(0oo2(2',O>  +  ---)  J 

=  e2  (01OO(2,  <)  0ioo(2',  f  )>  +  /2  (0oio(2, 0  0oio(2',  f  )>  + 

C2  (0ooi  (2,  t)  0ooi  (2',  0)  +  2Ce  (0ioo (2,  t)  0ooi  (2',  0) 

=  e2  P 9i00eioo  (2,    2',  if)  +  P  Pooo,  So,,,  (2,    2',  t')  + 

C2       e001  (2,  *;  2',  f )  +  2Ce  Peiooeool  (2,  t;  z',  0  (3-22) 


In  the  previous  equation,  no  correlation  has  been  assumed  between  either 
the  rainfall  and  the  saturated  hydraulic  conductivity,  or  the  water  uptake  and  the 
saturated  hydraulic  conductivity.  Since  the  components  Po100o100  and  Pe01090io  are 
identical  to  Pelo0lo  and  F£i01e01  that  have  been  determined  in  the  previous  Chapter, 
only  the  components  Pg001  q001  and  Pg100  g00l  will  be  derived  here  following  the  same 
procedure  as  that  described  in  the  previous  Chapter. 

Pre- multiplying  equation  (3.21)  by  a  soil- water  content  perturbation,  0ooi(2',t')> 
at  location  z'  and  time  t'  independent  of  z  and  t  and  taking  the  expected  value,  leads 
to  the  governing  equation  of  Pe001e001, 
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(3.23) 


The  previous  equation  requires  determining  jP0ooiX.  Post-multiplying  equation 
(3.21)  by  x(^),  at  location  z  and  time  t,  independent  of  z'  and  t'  and  taking  the 
expected  value,  leads  to  the  governing  equation  of  Pe00lX: 
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(3.24) 


where  Pxx  denotes  the  temporal  covariance  of  the  water  uptake  perturbation. 
The  governing  equation  of  Pe001e100  is  obtained  by  taking  the  expected  value  of  equa- 
tion (2.13)  derived  in  the  previous  Chapter  (and  that  is  equivalent  to  the  equation  of 
#ioo)  after  pre-multiplying  by  9ooi(z',t'),  at  a  location  z'  and  a  time  t',  independent 
of  z  and  t: 
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The  governing  equation  (3.25)  requires  determining  Pe0Q1s-  Therefore,  post- 
multiplying  equation  (3.21)  by  5(z  =  0,t),  at  location  z=0  and  time  t,  independent 
of  z'  and  t'  and  taking  the  expected  value,  leads  to  the  governing  equation  of  Pg001s'- 
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Darcy's  flux  statistics 

The  recharge  rate  from  the  unsaturated  zone  to  groundwater  is  controlled  by 
the  amount  of  water  infiltrating  the  soil  surface  as  well  as  by  the  amount  of  water 
uptake.  The  spatio-temporal  variability  of  the  recharge  rate  can  be  characterized  by 
the  Darcy's  flux  statistics.  Like  for  the  soil-water  content,  the  mean  and  the  variance 
of  the  Darcy's  flux  are  the  ensemble  moments  of  interest.  Inserting  the  power  series 
expansion  of  0,  and  the  Taylor's  series  expansion  of  T>(9)  and  K(9)  into  equation 
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(3.2)  for  the  Darcy's  flux  and  noting  that  ^ff-  =  0  leads  to  the  following  second 
order  equation: 
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(3.27) 


Taking  the  expected  value  of  equation  (3.27)  leads  to  the  following  equations 
for  the  zero,  first  and  second  order  estimates  of  the  mean  Darcy's  flux: 
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where  the  expressions  (^100^),  (floio^f),  and  (0Ooi  ^ff)  simplify  to 
^,  and  ^^ffS  respectively,  the  component  -P^ooo)^  +%\inM  ^010+^(^000) 


2     3z    >  2     9z     »  .-v,   "f— — "      "  vuuuy    ^     1  30I0QOO 

equals  to  zero  for  the  one-dimensional  flow  (equation  (2.18)  in  the  previous  Chapter) 
and  the  cross-covariances  Pv$  and  P^  are  assumed  to  be  zero.  Note  that  the  second 
order  mean  Darcy's  flux  can  be  calculated  once  &el00,  cr|010,  ^j001,  (#101),  (#200),  (#002), 
Pri9010(z  =  z'),  and  Pe100oool(z  =  z')  are  known.  The  second  order  dimensionless 
Darcy's  flux  autocovariance  can  be  approximated  as: 
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where  it  was  assumed  that  Pir]  =  0  and  P7)C  =  0.  It  is  worth  noting  that  i) 
because  of  the  specification  of  the  bottom  boundary  condition  (equation  (2.18)  in 
the  previous  Chapter)  the  autocovariance  of  the  Darcy's  flux  depends  only  upon  the 
variabilities  of  the  rainfall  and  water  uptake,  and  ii)  the  autocovariance  of  the  Darcy's 
flux  includes  equation  (2.47)  and  two  new  terms  which  are  functions  of  (2  and  (e  that 
capture  the  effect  of  the  water  uptake. 

Pore-water  velocity  statistics 

Rainfall  and  water  uptake  control  the  amount  of  solute  discharge  to  groundwa- 
ter as  well  as  the  travel  time  of  the  solute  in  the  vadose  zone.  Since  infiltrating  water 
carries  inert  and  reactive  solute  through  soils,  it  is  important  to  quantify  the  spatio- 
temporal  variability  of  the  pore- water  velocity  due  to  the  spatially  random  structure 
of  the  saturated  hydraulic  conductivity,  the  temporally  random  structure  of  rainfall 
series  and  the  spatially  variable  and  temporally  distributed  random  structure  of  water 
uptake.  Our  interest  focuses  on  the  mean  and  the  variance  of  the  pore-water  velocity, 
v(z,t)  defined  as  the  ratio  of  two  random  processes, 


q(9(z,t)) 

In  order  to  derive  the  statistics  of  v,  an  expansion  of  the  previous  equation 
about  (?ooo  is  required  {Mood  et  ai,  1974), 
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Substituing  equations  (3.6)  and  (3.27)  into  equation  (3.29)  leads  to 
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(3.30) 


where  it  was  assumed  that  Pvs  =  0  and  P^x  —  0.  Taking  the  expected  value  of 
equation  (3.30)  produces  the  following  zero,  first  and  second  order  estimates  of  the 
mean  pore-water  velocity: 
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where  the  same  simplifications  and  assumptions  used  to  derive  the  mean 
Darcy's  flux  have  been  applied.  The  second  order  dimensionless  pore-water  veloc- 
ity autocovariance  can  be  approximated: 
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(3.31) 


where  it  was  assumed  that  Psv  —  P^  =  0.  Compared  to  the  Darcy's  flux 
autocovariance,  it  is  worth  noting  that  i)  the  autocovariance  of  the  pore- water  velocity 
depends  upon  the  variabilities  of  the  rainfall,  of  the  water  uptake,  as  well  as  that  of 
the  saturated  hydraulic  conductivity  and  ii)  the  pore-water  velocity  autocovariance 
includes  equation  (2.51)  and  two  new  terms  which  are  functions  of  C2  and  that 
capture  the  effect  of  water  uptake. 
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3.2.3    Semi-analytical  Solutions  of  the  Flow  Ensemble  Moments 

All  the  ensemble  moments,  Pqq,  Pvv,  Pee  are  a  function  of  P01OO01OO,  ^oiofloio 
PflooiOoon  and  Fe001e100.  The  derivations  and  solutions  of  Pe100e100,  ^Woio  have  been 
completed  in  the  previous  Chapter.  Therefore,  only  the  expressions  for  Pe001e00i  and 
Peooie100  must  be  derived.  To  solve  the  equation  for  P9ooio001,  a  double  Laplace  trans- 
form is  applied  over  t  and  t'  in  equations  (3.23)  and  its  respective  boundary  and 
initial  conditions, 
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where  the  double  tilde  represents  the  double  Laplace  transform  and  s  and  s' 
are  the  Laplace  variables  corresponding  to  t  and  t',  respectively.  The  solution  of  this 
second  order  linear  and  nonhomogeneous  differential  equation  is  given  by: 
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and  ^  =  /l  +  Ag^l6oo°l2 .  The  previous  solution  requires  the  solution  of  Pe00lX. 
Again,  a  double  Laplace  transform  is  applied  over  t  and  t'  in  equations  (3.24)  and  its 
respective  boundary  and  initial  conditions 
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The  solution  of  this  second  order  linear  and  nonhomogeneous  differential  equa- 
tion is  given  by: 
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and  ^2  =  Jl+  4^y°°o)2  •  Similarly,  taking  the  double  Laplace  transform  of 
equation  (3.25)  and  its  respective  boundary  and  initial  conditions  gives: 
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The  solution  of  this  second  order  linear  and  homogeneous  differential  equation 
is  given  by: 
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and  &  is  denned  as  before.  Note  that  in  order  to  derive 
required.  Taking  the  double  Laplace  Transform  of  equation  (3.26)  and  its  boundary 
and  initial  conditions  yields: 
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The  solution  of  this  second  order  linear  and  nonhomogeneous  differential  equa- 
tion is  given  by: 
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and  £c  is  defined  as  before.  Finally,  the  components  of  the  soil-water  content 
autocovariance  derived  in  the  previous  Chapter,  Pe100ei0o  and  ^0oio0oio>  35  we^  as 
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P0OO10OO1  and  Pftufao  can  then  be  substituted  into  the  double  Laplace  transform  of 
equation  (3.22): 

Pio{z\  s'\  z,  s)  =  e2  P^ih00(z',  s';  z,  s)  +  f2  Po01080io(z' >  s'>  z> s)  + 

c2  K^L  (2'> s';  *. s)  +  2^e  *>  ^)  (3.40) 

The  autocovariances  of  the  Darcy's  flux  and  pore-water  velocity  are  obtained 
by  taking  the  double  Laplace  transform  of  equations  (3.28)  and  (3.31): 
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A  numerical  scheme  was  implemented  to  invert  the  Laplace  transforms  Pee, 
P^q,  and  P»v  into  real  time.  Taking  advantage  of  the  relationship  between  the  Laplace 
and  the  Fourier  transforms  (Appendix  1),  the  two-dimensional  Fast  Fourier  Transform 
algorithm  of  Press  et  al,  (1992)  was  used  to  obtain  Pe6(z' ,  t';  z,  t),  Pqq{z',  If  ;  z,  t),  and 
Pvv(z',t';z,t).  The  dimensional  autocovariances  of  the  Darcy's  flux  and  pore-water 
velocity  in  real  time  are  obtained  by  inverting  equations  (3.41)  and  (3.42)  using 
the  two-dimensional  Fast  Fourier  Transform  algorithm  and  then  multiplying  by  the 
scaling  factor  ICq2. 
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3.3    Illustration  of  Theoretical  Results 

The  means  and  the  autocovariances  of  the  soil- water  content,  Darcy's  flux  and 
pore- water  velocity  are  illustrated  for  a  heterogeneous  fine  sandy  soil  subject  to  ran- 
dom series  of  rainfall  and  water  uptake  representative  of  Tampa,  Florida.  Under  such 
climatic  conditions,  the  mean  volume  of  water  removed  by  water  uptake  represents  60 
%  of  the  mean  rainfall.  In  the  example,  the  characteristic  rooting  depth  is  assumed  to 
be  1  meter  (z*  =  1  m).  The  soil-moisture  diffusivity  and  relative  hydraulic  conductiv- 
ity are  modeled  using  the  Brooks  and  Corey  relationships  (Brooks  and  Corey,  1964)  ■ 
This  empirical  model  is  well  adapted  for  sandy  soils  (Lappala  el  ai,  1993 )  and  usually 
fits  measured  field  data  better  than  the  less  complex  Gardner-Russo  model  (Zhang 
et  ai,  1998).  The  Brooks-Corey  model  represents  the  soil-moisture  diffusivity  and 
relative  hydraulic  conductivity  as  follows: 

where  9S  is  the  saturated  soil-water  content,  assumed  equivalent  to  the  poros- 
ity, 9r  the  residual  soil-water  content,  and  A  the  grain  size  distribution  parameter. 
The  Taylor's  series  expansion  of  V  and  K  about  0OOo  are  expressed  as: 
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Table  3.1  shows  the  ensemble  moments  assumed  for  the  rainfall  and  water 
uptake  series  and  for  the  saturated  hydraulic  conductivity  while  Table  3.2  presents 
the  deterministic  properties  of  the  soil. 

In  this  example,  the  variance  of  the  saturated  hydraulic  conductivity  was  se- 
lected as  the  upper  bound  of  /  so  that  it  would  allow  a  converging  power  series 
solution  for  9.  The  statistics  of  the  rainfall  were  calculated  from  long  term  weather 
records  from  Tampa,  Florida,  and  the  mean  water  uptake  was  set  to  be  equal  to  60 
%  of  the  mean  rainfall. 

3.3.1    Mean  of  the  Soil- water  Content.  Darcv's  Flux  and  Pore- water  Velocity 

After  specifying  the  hydraulic  conductivity  relationship,  the  equivalent  steady- 
state  soil- water  content  is  determined  by  solving  equation  (3.20): 


R  - £  -  <TSr>i+3  = 0 

V  s  —  V  r 


This  nonlinear  equation  is  solved  using  the  Newton-Ralphson  method.  Then, 
the  first  order  dimensional  means  of  the  Darcy's  flux  and  of  the  pore-water  velocity 
are  found  to  be  equal  to: 
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where  the  dimensionless  quantities  q  and  v  are  multiplied  by  the  scaling  factor 
JC*G.  The  power  series  of  9  converges  rapidly  due  to  the  small  values  of  e  =  1.6  x  10-5, 
/  =  1  and  (  =  6  x  10"7.  Thus,  the  calculations  of  the  second  order  correction 
to  the  mean  soil- water  content,  Darcy's  flux  and  pore- water  velocity  have  not  been 
undertaken. 
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3.3.2    Autocovariances  of  the  Soil-water  Content.  Darcy's  Flux  and  Pore-water  Velocity 

Evaluating  the  autocovariances  of  the  soil-water  content,  Darcy's  flux  and 
pore-water  velocity  requires  specifying  temporal  correlation  structures  for  the  rainfall, 
and  for  the  water  uptake,  as  well  as  a  spatial  correlation  structure  for  the  saturated 
hydraulic  conductivity.  Here  the  auto-covariance  for  water  uptake  and  the  cross- 
covariance  between  the  water  uptake  and  the  rainfall,  P6x,  were  derived  from  assuming 
that  the  water  uptake  is  a  weighted  combination  of  two  fields,  one  correlated  to  the 
rainfall  and  one  uncorrelated  (Gelhar,  1993),  that  is 


X  =  Psx6  +  \jl-Psxv 
where  v  is  a  random  zero-mean  temporally  correlated  process  that  is  uncorre- 
lated with  5,  and  pSx  is  a  weighting  parameter  (-1  <  psx  <  1)  that  determines  the 
correlation  between  x  and     Tne  covariances  used  in  this  Chapter  are  presented  in 
Table  3.3. 

The  autocovariances  of  the  soil-water  content,  Darcy's  flux  and  pore-water 
velocity  are  evaluated  at  two  levels  of  cross-correlation  (psx  —  1,  and  -  1)  that  are 
not  necessarily  realistic  but  bound  possibilities  of  correlation  and  the  resulting  re- 
sponse of  the  flow  system.  However,  the  positively  correlated  case  can  be  interpreted 
in  the  context  of  an  active  vegetation  removing  increased  volumes  of  water  around 
rainfall  events.  A  negatively  correlated  case  can  be  interpreted  as  modeling  reduced 
evapotranspiration  due  to  increased  cloud  cover  near  precipitation  events. 

3.3.3    Discussion  of  Theoretical  Results 

The  equivalent  steady-state  mean  soil- water  content  and  the  coefficients  of  the 
soil-moisture  diffusivity  and  relative  hydraulic  conductivity  functions  are  presented 
in  Table  3.4.  The  anticipated  effect  of  the  water  uptake  is  confirmed.  Away  from 
the  land  surface,  the  equivalent  steady-state  mean  soil-water  content  is  smaller  under 
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the  water  uptake  scenario  than  with  no  uptake.  This  implies  that  the  mean  Darcy's 
flux  and  the  mean  pore-water  velocity  are  also  reduced.  Note  that  the  true  mean 
soil- water  content  decreases  with  depth  from  .100858  at  the  surface  to  the  equivalent 
steady-state  mean  value  of  .092  deep  in  the  profile.  The  mean  soil-water  content  is 
reached  at  a  depth  controlled  by  the  magnitude  of  the  characteristic  rooting  depth 
z*. 

The  effect  of  the  water  uptake  over  depth  on  the  second  order  variances  of 
soil-water  content,  Darcy's  flux  and  pore-water  velocity  in  a  fine  sandy  soil  are  shown 
in  Figure  3.1.  The  first  case  (Case  1)  does  not  include  any  water  uptake  (and  thus  it 
is  equivalent  to  the  case  examined  in  Chapter  2),  the  second  case  (Case  2)  considers 
a  deterministic  water  uptake  (a\.  =  0)  equivalent  to  the  mean  water  uptake  of  the 
two  last  cases  (Case  3  and  4)  where  60  %  of  the  mean  infiltrating  water  is  removed  by 
the  mean  water  uptake  at  two  different  levels  of  cross-correlation.  This  figure  shows 
that  i)  at  depths  greater  than  the  characteristic  rooting  depth  (where  the  equivalent 
steady-state  mean  soil- water  content  is  approached),  the  effect  of  the  water  uptake  is 
to  increase  the  variance  of  the  soil- water  content  and  to  reduce  that  of  the  Darcy's  flux 
and  pore-water  velocity  compared  to  no  water  uptake,  ii)  the  statistics  predicted  with 
a  deterministic  water  uptake  are  positioned  between  those  of  the  two  levels  of  cross- 
correlation,  iii)  the  variances  of  the  soil-water  content,  Darcy's  flux  and  pore-water 
velocity  undergo  a  sharper  decrease  over  depth  with  water  uptake  than  without,  and 
iv)  compared  to  a  positive  correlation  between  rainfall  and  water  uptake,  a  negative 
p5x  increases  the  variances  of  soil- water  content,  Darcy's  flux  and  pore- water  velocity. 
Results  from  a  previous  investigation  assuming  steady  state  flow  by  Rubin  and  Or 
(1993)  agree  with  these  results  observed  for  the  soil-water  content.   Among  their 
investigations,  Rubin  and  Or  examined  the  effect  of  an  exponential  water  uptake 
function  with  constant  transpiration  rate  on  the  spatial  evolution  of  pressure  variance 
in  a  vertical  bounded  domain.  They  observed  that  an  increase  of  the  water  uptake 
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depth  or  an  increase  of  the  transpiration  leads  to  a  pressure  variance  that  is  larger 
than  that  observed  without  water  uptake,  and  that  the  water  uptake  cause  steeper 
negative  pressure  gradients  near  the  soil  surface. 

The  effect  of  the  water  uptake  on  the  pore-water  velocity  has  been  previously 
investigated  by  Russo  et  al.  (1998).  They  performed  a  three-dimensional  single 
numerical  simulation  of  unsaturated  flow  in  heterogeneous  conditions  with  a  non- 
monotonic but  deterministic  rainfall  pattern.  Their  water  uptake  was  expressed  as  a 
function  of  the  local  pressure  and  hydraulic  conductivity  and  was  applied  throughout 
the  domain.  The  authors  found  that  the  water  uptake  suppresses  high  velocity  and 
consequently  reduces  the  velocity  variance  in  the  longitudinal  direction.  They  also 
showed  that  the  soil-water  content  variance  is  enhanced  by  water  uptake.  Such  nu- 
merical results  support  and  agree  with  the  theoretical  development  and  predictions 
of  the  current  stochastic  flow  model.  So  far,  the  effect  of  water  uptake  on  the  Darcy's 
flux  variance  has  not  yet  been  reported.  However,  since  the  evolution  of  its  variance 
over  depth  behaves  similarly  to  that  of  the  pore-water  velocity,  it  is  most  likely  that 
the  water  uptake  also  suppresses  high  values  of  Darcy's  flux  and  therefore  reduces  its 
variability. 

The  effect  of  the  cross-correlation  between  rainfall  and  water  uptake  on  the 
variances  of  soil-water  content,  Darcy's  flux  and  pore-water  velocity  can  be  explained 
as  either  reducing  or  enhancing  the  effect  of  rainfall  on  the  same  statistics.  From 
the  governing  equation  (3.1),  the  rainfall  and  the  water  uptake  have  opposite  effects 
on  the  water  budget  in  the  soil,  the  former  adds  water  into  the  soil  while  the  latter 
removes  it.  This  relationship  is  actually  controlled  by  the  sign  of  the  cross-correlation. 
A  positive  pSx  means  that  large  rainfall  events  coincide  with  large  water  uptake  values 
that  therefore  remove  water  from  the  soil.  On  the  contraty,  a  negative  psx  means  that 
large  rainfall  events  coincide  with  small  water  uptake  values,  therefore  increases  the 
water  in  the  soil.  Therefore,  the  larger  variances  of  soil-water  content,  Darcy's  flux 
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and  pore-water  velocity  observed  in  the  case  of  a  negative  cross-correlation  are  due  to 
the  effect  of  the  rainfall  and  water  uptake  variabilities  working  together  to  enhance 
water  content  variability,  whereas  in  the  case  of  a  positive  cross-correlation,  rainfall 
and  water  uptake  work  against  each  other  to  reduce  water  content  variability. 

Figures  3.2-3.3  depict  the  normalized  temporal  autocovariances  of  the  soil- 
water  content  in  a  heterogeneous  fine  sandy  soil  at  two  different  depths.  The  autoco- 
variances were  normalized  by  dividing  by  the  variance  calculated  at  the  soil  surface. 
Like  in  the  previous  Chapter,  these  figures  show  asymptotic  temporal  stationarity  of 
the  soil- water  content  autocovariance  which  is  reached  rapidly  near  the  surface.  These 
figures  show  also  that  soil-water  content  is  correlated  over  greater  times  as  depth  in- 
creases. The  stationary  functions  used  to  model  the  covariances  of  the  rainfall  and 
water  uptake  series  explain  the  stationary  behavior  of  the  soil-water  content  over 
time.  The  propagation  of  the  soil- water  content  through  soil  subject  to  water  uptake 
is  governed  by  a  set  of  nonhomogemeous  partial  differential  equations  that  resemble 
advection-dispersion  equations  with  an  internal  sink  term.  In  each  of  these  equations, 
the  soil-moisture  diffusivity  produces  a  diffusive  effect  in  the  soil  and  explains  why 
the  soil-water  content  has  smaller  variance  and  is  more  correlated  over  time  as  depth 
increases.  Figures  3.4-3.5  present  the  normalized  temporal  autocovariances  of  the 
Darcy's  flux  and  pore-water  velocity  in  a  heterogeneous  fine  sandy  soil  at  0.5  meter. 
Since  the  normalized  autocovariances  of  the  Darcy's  flux  and  pore-water  velocity  are 
controlled  by  the  autocovariance  of  the  soil- water  content,  their  temporal  behavior 
is  similar  to  that  of  the  soil-water  content  and  is  therefore  only  demonstrated  at  one 
depth. 

Figures  3.6-3.7  display  the  normalized  spatial  autocovariances  of  the  soil-water 
content  in  a  heterogeneous  fine  sandy  soil  at  two  different  observation  times.  These 
figures  confirm  the  temporal  stationarity  of  the  soil-water  content  and  indicate  its 
non-stationary  evolution  over  space.  Once  again,  the  soil-moisture  diffusivity  in  the 
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flow  governing  equations  explains  the  non-stationarity  of  the  soil-water  content  over 
space.  Figures  3.8-3.9  present  the  normalized  spatial  autocovariances  of  the  Darcy's 
flux  and  pore-water  velocity  in  a  heterogeneous  fine  sandy  soil  at  0.5  meter.  These 
surface  plots  show  also  the  non-stationary  behavior  over  depth  of  the  Darcy's  flux 
and  pore-water  velocity. 

Taking  advantage  of  the  additive  properties  of  the  effect  of  each  random  pro- 
cess in  the  statistics  of  soil-water  content,  Darcy's  flux  and  pore-water  velocity,  it 
is  possible  to  evaluate  the  relative  contribution  of  each  of  them  on  the  overall  un- 
certainty. This  analysis  is  performed  on  the  pore-water  velocity  variance  for  the 
four  cases  mentioned  earlier  at  two  different  depths  and  is  summarized  in  Table  3.5. 
These  results  show  i)  a  dominance  of  the  rainfall  variability  over  those  of  the  satu- 
rated conductivity  and  water  uptake  at  any  depth,  ii)  the  increasing  contribution  of 
the  saturated  conductivity  over  depth,  iii)  the  decreasing  contribution  of  water  up- 
take over  depth,  and  iv)  a  reduction  or  an  increase  of  the  water  uptake  contribution 
depending  upon  the  sign  of  pSx.  Note  that  i)  the  contribution  of  the  saturated  hy- 
draulic conductivity  is  larger  without  water  uptake  because  the  constant  coefficients 
{V(§00o)  and  %\-9oJ  of  the  governing  flow  equation  are  larger,  and  ii)  that  below  the 
characteristic  water  uptake  depth  the  contribution  of  the  water  uptake  is  zero.  Since 
the  equivalent  steady-state  mean  soil- water  content  is  not  valid  at  .1  meter,  the  cal- 
culations made  at  that  depth  are  not  as  reliable  as  those  made  at  5  meters.  However, 
the  relative  magnitude  of  these  quantities  help  to  explain  the  relative  contribution  of 
each  random  process  on  the  overall  variability. 

So  far,  we  have  noted  that  the  variability  in  pore-water  velocity  is  dominated 
by  the  contribution  of  the  rainfall  variability  that  decays  over  depth.  Without  loss  of 
generalization,  the  same  observation  can  be  established  also  for  the  Darcy's  flux  and 
the  soil- water  content  variances.  Unlike  the  soil- water  content,  we  have  also  noted  that 
the  effect  of  the  water  uptake  is  to  reduce  the  variance  of  the  pore- water  velocity  and 
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Darcy's  flux  in  the  region  where  the  assumption  of  an  equivalent  steady-state  mean 
soil-water  content  holds.  The  next  analysis  will  focus  on  the  Darcy's  flux  autocovari- 
ance  to  shed  some  light  on  why  the  Darcy's  flux  variance  is  smaller  with  water  uptake 
compared  to  no  water  uptake.  Since  the  rainfall  is  the  dominant  source  of  uncertainty, 
only  the  components  of  e2  in  equation  (3.28)  are  going  to  be  investigated.  The  com- 
ponents of  the  Darcy's  flux  autocovariance  are  ( f^J^Wioo  defined  as  the  gravi- 
tational component,  and  (^ooo))2^^  and  V(9000)  f  b000  + 
defined  as  the  two  dispersive  components.  Each  component  has  been  calculated  for 
the  Darcy's  flux  variance  at  two  different  depths  with  no  water  uptake  (Case  1)  and 
with  water  uptake  positively  correlated  with  rainfall  (Case  4).  The  results  are  sum- 
marized in  Table  3.6. 

For  any  given  case,  the  Darcy's  flux  autocovariance  experiences  a  redistribution 
over  depth  of  its  components.  At  any  given  depth  the  gravitational  component  largely 
exceeds  the  two  dispersive  components,  and  more  particularly  at  deeper  depths.  For 
the  two  cases,  the  contribution  of  the  two  dispersive  components  is  reduced  over 
depth  while  that  of  the  gravitational  component  is  increased.  Moreover,  at  shallow 
depths,  the  effects  of  water  uptake  are  to  enhance  the  contribution  of  the  dispersive 
components  and  to  reduce  the  contribution  of  the  gravitational  component. 

The  dispersive  components  are  the  product  of  the  spatial  derivatives  of  the 
soil- water  content  variance  multiplied  by  the  constants  %\qqoo  and  £>(0Ooo)  which 
are  evaluated  at  the  equivalent  steady-state  mean  soil- water  content.  At  the  shallow 
depth  of  0.1  meter,  the  sharper  gradient  of  the  soil- water  content  variance  due  to 
water  uptake  causes  an  increase  of  the  spatial  derivatives  +  ^igftm )  and 

d2Pdldz'm  tnat  is  !arger  than  the  reduction  of  the  multiplicative  constants  which  are 
evaluated  at  a  smaller  0Ooo  compared  to  no  water  uptake  (Table  3.4).  As  a  result,  the 
water  uptake  increases  the  dispersive  components  of  the  Darcy's  flux  autocovariance. 
Concerning  the  gravitational  component,  the  effect  of  water  uptake  is  to  reduce  %\gfm 
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more  than  P9l0Qe100  is  increased.  As  a  result,  the  gravitational  component  is  reduced 
because  of  the  water  uptake.  The  same  explanations  apply  also  at  the  depth  of  1.5 
meter.  Therefore,  the  effect  of  the  water  uptake  is  a  decrease  of  the  gravitational 
component  of  the  Darcy's  flux  autocovariance  that  is  larger  than  the  increase  of  the 
dispersive  components  so  that  a  smaller  Darcy's  flux  variance  is  produced. 

The  effects  of  the  water  uptake  on  temporal  and  spatial  autocorrelation  of  the 
soil- water  content,  Darcy's  flux  and  pore- water  velocity  are  shown  on  Figures  3.10- 
3.13  for  the  two  cross-correlation  levels  between  rainfall  and  water  uptake  and  two 
different  depths.  The  time  and  depth  of  observation  have  been  chosen  so  that  i)  the 
initial  time  occurs  when  temporal  stationarity  is  satisfied  and  ii)  the  assumption  of 
equivalent  steady-state  mean  soil-water  content  is  valid  at  the  depth  of  observation. 
The  effect  of  the  water  uptake  at  a  given  depth  is  shown  to  increase  the  temporal 
correlation  scale  of  the  soil-water  content,  Darcy's  flux  and  pore-water  velocity  (Fig- 
ures 3.10  and  3.11),  but  to  have  a  limited  effect  on  the  spatial  correlation  (Figures 
3.12  and  3.13).  Moreover,  the  positive  cross-correlation  level  slightly  enhances  this 
effect  compared  to  the  negative  one.  Both  the  temporal  and  spatial  correlation  scales 
increase  with  depth. 

3.4  Summary 

In  this  Chapter,  a  one-dimensional,  stochastic  model  was  developed  to  describe 
water  flow  in  a  variably  saturated  and  heterogeneous  porous  medium  subject  to  a 
transient  flow  regime  resulting  from  temporally  random  rainfall  and  spatially  variable 
and  temporally  random  water  uptake.  The  model  predicts  the  ensemble  moments  of 
the  soil-water  content,  Darcy's  flux,  and  pore- water  velocity  in  a  semi-analytical  form 
as  a  function  of  the  statistical  parameters  of  the  rainfall,  water  uptake,  and  saturated 
hydraulic  conductivity.  The  one-dimensional  soil- water  content  mean  and  covariances 
equations  are  developed  by  assuming  a  converging  series  solution  for  the  soil-water 
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content.  Approximate  ensemble  flow  moment  equations  are  derived  in  Laplace  space 
and  are  inverted  using  a  double  Fast  Fourier  Transform  code.  The  predictions  of  the 
developed  model  were  illustrated  for  a  heterogeneous  fine  sandy  soil  subject  to  rainfall 
and  water  uptake  series  recorded  in  Florida,  and  were  compared  to  previous  works. 

The  main  findings  concerning  the  effect  of  the  water  uptake  on  the  water  flow 
statistics  are  summarized  as  follows: 

-  Compared  to  a  no  uptake  scenario,  water  uptake  increases  the  spatial  gradient 
of  the  soil-water  content  variance  at  shallow  depths.  Moreover,  the  soil-water  content 
variances  are  enhanced  at  all  depths  by  both  the  random  and  deterministic  water 
uptake,  regardless  of  the  cross-correlation  assumption  between  rainfall  and  water 
uptake.  Consideration  of  random  water  uptake  reduces  the  variability  of  the  Darcy's 
flux  and  pore-water  velocity  at  all  depths  compared  to  a  no  water  uptake  scenario 
for  any  cross-correlation  assumption.  However,  a  negative  cross-correlation  case  has 
higher  Darcy's  flux  and  pore-water  velocity  variance  over  depth  compared  to  the 
positive  cross-correlation  case,  and  appears  to  converge  to  the  no  uptake  case  at 
deeper  depths  for  the  pore- water  velocity  variance. 

-  The  effect  of  the  water  uptake  is  shown  to  increase  the  temporal  correlation 
scales  of  the  soil- water  content,  Darcy's  flux  and  pore- water  velocity  for  all  cross- 
correlation  assumptions  but  to  have  a  limited  impact  on  the  spatial  correlation  of 
these  variables.  The  effect  on  temporal  correlation  is  amplified  over  depth  due  to  the 
advection-dispersion  form  of  the  flow  equation.  It  is  worth  noting  that  the  soil- water 
content  autocovariance  controls  that  of  the  Darcy's  flux  and  pore-water  velocity  and 
and  that  the  rainfall  variability  dominates  that  of  the  water  uptake  and  saturated 
hydraulic  conductivity  in  determining  the  overall  uncertainty. 

-  Before  comparing  the  current  model  to  numerical  simulations  or  field  data, 
some  improvements  are  required.  First,  in  the  current  work,  the  non-linear  decay  of 
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the  mean  soil-water  content  over  depth  is  replaced  by  a  constant  equivalent  steady- 
state  value  that  is  valid  only  in  the  region  of  the  vadose  zone  where  the  effect  of  the 
time- varying  boundary  conditions  are  strongly  damped.  Relaxing  this  assumption  will 
lead  to  partial  differential  flow  equations  with  nonconstant  coefficients  that  must  be 
solved  using  numerical  techniques.  Second,  the  current  model  would  be  more  realistic 
if  the  water  uptake  were  a  function  of  the  local  soil- water  content.  The  solution 
technique  required  to  solve  this  new  equation  will  depend  upon  the  complexity  of 
the  relationship.  A  water  uptake  flux  expressed  as  a  linear  function  of  the  local  soil 
water  could  be  solved  with  the  current  methodology  whereas  a  non-linear  function 
may  require  a  numerical  solution.  However,  the  dissipating  variability  of  the  soil- 
water  content  over  depth  suggests  that  the  effect  of  a  soil-water  content-dependent 
water  uptake  flux  on  the  statistics  of  soil-water  content,  Darcy's  flux  and  pore-water 
velocity  would  resemble  the  results  of  the  current  model. 
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Table  3.1:  Statistical  properties  of  the  rainfall,  the  water  uptake  and  the  log-saturated 
hydraulic  conductivity. 


Variable 

Mean 

Variance 

Correlation  length 

Rainfall 
Ln  /C; 

Water  uptake 

(q*)  =  0.0038  m/d 

K,q  —  5  m/d 
(e*)  =  0.0023  m/d 

al  =  9.81  x  10"5m2/d2 
afnK,  =  1 
al  =  4  xlO-W/d2 

Xs  =  .6  d 
A,,  =1  m 
Ax  =.6  d 

Table  3.2:  Deterministic  hydraulic  properties  of  the  fine  sandy  soil. 


Property 

Value 

.35 

Or 

.063 

K 

-.82  m 

X 

3.7 

Table  3.3:  Covariances  used  in  this  Chapter  and  their  Laplace  or  Fourier  transforms. 


Covariance 

Expressions 

Laplace/Fourier  transform 

Pss 
p 

p 
p 

|*— *'| 

e  x« 

P]  P55  +  (1  -  p2sx)P»v 

e  Xji 

n5[S,b)  -  (i+xss'){s+s')(l+\ss) 

i   -\  K(2+K(s+s')) 
rvv\a  ,!>)  —  (i+A^s'Xs+s'Xl+A^) 

tj>X}S  '  S>  ~  (l+A4y)(s+s')(l+*<ss)   

P^(s'-,s)  =  p2sxP^s  +  (l-p2sx)P7u 

2  \ 

DVV  -  7r(l+A2fcJ) 

115 


Table  3.4:  Equivalent  steady-state  mean  soil- water  content,  and  its  associated  zero 
order  soil-moisture  diffusivity  and  first  derivative  of  the  hydraulic  conductivity. 


Quantity 

No  water  uptake 

60  %  water  uptake 

#000 

.100858 

.0920 

dK  1 

.359073 

.1822 

d6  Uooo 

T>(6000) 

.0388 

.0212 

<7*(0ooo)  (m/d) 

.0038 

.0015 

v*(0m){m/d) 

.0380 

.0162 

Table  3.5:  Effect  of  the  water  uptake  at  two  different  depths  on  the  magnitude  of  each 
component  of  the  pore- water  velocity  variance  (z=z\  in  meters).  The  percentages 
denote  the  contribution  of  each  component  to  the  pore-water  velocity  variance. 


Case  1 

Case  2 

Case  3 

Case  4 

(xlO-4 
m2/d2) 

No  water 
uptake 

Constant 
60%,  a\.  =  0 

Random 
60%,pSx  =  -1 

Random 
60%,  pSx  =  1 

.1  m 

5.  m 

.1  m 

5.  m 

.1  m 

5.  m 

.1  m 

5.  m 

2  p 

t  -'01000100 

65.88 

(99.8%) 

11.59 

(99%) 

68.91 

(99.97%) 

6.88 

(99.71%) 

68.91 

(93.34%) 

6.88 

(76.2%) 

68.91 

(107.33%) 

6.88 

(144.23%) 

2^-^0010100 

0 

(0%) 

0 

(0%) 

0 

(0%) 

0 

(0%) 

4.81 

(6.52%) 

2.13 

(23.6%) 

-4.81 

(-7.5%) 

-2.13 

(-44.65%) 

C  -ffoO 1^001 

0 

(0%) 

0 

(0%) 

0 

(0%) 

0 

(0%) 

.08 

(.11%) 

0. 

(0%) 

.08 

(.13%) 

0. 

(0%) 

^-Pflooifooi 

.14 

(.2%) 

.14 

(i%) 

.02 

(.03%) 

.02 

(.29%) 

.02 

(.03%) 

.02 

(.2%) 

.02 

(.03%) 

.02 

(.42%) 

9 

67.35 

11.73 

68.90 

6.90 

73.82 

9.03 

64.2 

4.77 
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Table  3.6:  Effect  of  the  water  uptake  at  two  different  depths  on  the  magnitude  of  each 
component  of  the  Darcy's  flux  variance  (z=z\  in  meters).  The  percentages  denote 
the  contribution  of  each  component  to  the  Darcy's  flux  variance. 


Case  1 

Case  4 

Quantity 

No  water 
uptake 

60% 
£0=.092,  pSx  =  1 

(xl0-6m2/d2) 

z  =  .1  m 

z  =  5  m 

z  =.1  m 

z  =  5  m 

(2>(4o))2^^ 

6.48 
(8%) 

.084 

(.6%) 

9.7 

(14%) 

.047 
(.6%) 

P(0ooo)  Hb000 

/■ap«ioo«ioo    |   Psioo«ioo  \ 
I       dz         1        dz'  I 

9.73 
(12%) 

.176 
(1.2%) 

12.32 
(18%) 

.10 

(1.5%) 

(dK\  \2pn 

64.80 
(80%) 

14.4 
(98.2%) 

45.13 
(67%) 

6.84 
(97.7%) 

Sum 

81.01 

14.66 

67.08 

7.04 
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c)  Pore  water  velocity 


Depth,  in  meters 


 60%  water  uptake,  (^^=1 

 60%  water  uptake,  f^^=-l 

—  -  -  60%  constant  water  uptake 

—  —  No  water  uptake 


Figure  3.1:  Effect  of  the  water  uptake  on  the  variance  over  depth  of  a)  the  soil- water 
content,  b)  the  Darcy's  flux  and  c)  the  pore-water  velocity. 
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0.0  0.2  0.4  0.6  0.8  1.0 
Normalized  soil-water  content  autocovariance 


Figure  3.2:  Temporal  evolution  of  the  normalized  autocovariance  of  the  soil-water 
content  in  a  heterogeneous  fine  sandy  soil  with  60  %  mean  water  uptake  scenario  and 
pSx  =  1  at  z=z'=0.5  meter,  ^Zf^'s^  ■ 
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0.0  0.2  0.4  0.6  0.8  1.0 
Normalized  soil-water  content  autocovariance 


Figure  3.3:  Temporal  evolution  of  the  normalized  autocovariance  of  the  soil-water 
content  in  a  heterogeneous  fine  sandy  soil  with  60  %  mean  water  uptake  scenario  and 
pSx  =  1  at  z=z'=4  meters,  PM*=^,*;*'=W) 
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0.0       0.2       0.4        0.6       0.8  1.0 
Normalized  Dairy's  flux  autocovariance 


Figure  3.4:  Temporal  evolution  of  the  normalized  autocovariance  of  the  Darcy's  flux 
in  a  heterogeneous  fine  sandy  soil  with  60  %  mean  water  uptake  scenario  and  psx  —  1 
at  z=z'=.5  meter,  ^=f^'^m,t') 
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0.0       0.2        0.4       0.6       0.8  1.0 
Normalized  pore-water  velocity  autocovariance 


Figure  3.5:  Temporal  evolution  of  the  normalized  autocovariance  of  the  pore-water 
velocity  in  a  heterogeneous  fine  sandy  soil  with  60  %  mean  water  uptake  scenario  and 
pSx  =  1  at  z=z'=.5  meter,  ^J^fj^  ■ 


122 


z,  meters 


0.0  0.2  0.4  0.6  0.8  1.0 
Normalized  soil-water  content  autocovariance 


Figure  3.6:  Spatial  evolution  of  the  normalized  autocovariance  of  the  soil-water  con- 
tent in  a  heterogeneous  fine  sandy  soil  with  60  %  mean  water  uptake  scenario  and 
p6x  =  1  at  60  days,  P^^^ . 
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0.0  0.2  0.4  0.6  0.8  1.0 
Normalized  soil-water  content  autocovariance 


Figure  3.7:  Spatial  evolution  of  the  normalized  autocovariance  of  the  soil- water  con- 
tent in  a  heterogeneous  fine  sandy  soil  with  60  %  mean  water  uptake  scenario  and 
pSx  =  1  at  100  days,  ^'Z^M>. 
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2.0  z',  meters 


i     1     i     •     i     *     I  i     *  I 

0.0       0.2        0.4       0.6       0.8  1.0 

Normalized  Darcy's  flux  autocovariance 


Figure  3.8:  Spatial  evolution  of  the  normalized  autocovariance  of  the  Darcy's  flux  in 
a  heterogeneous  fine  sandy  soil  with  60  %  mean  water  uptake  scenario  and  psx  =  1 
at  60  days,  P"{%™*Z^°d)  ■ 
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Figure  3.9:  Spatial  evolution  of  the  normalized  autocovariance  of  the  pore-water 
velocity  in  a  heterogeneous  fine  sandy  soil  with  60  %  mean  water  uptake  scenario  and 
pSx  =  1  at  60  days,  P"l'^f^ . 
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a)  Soil  water  content  b)  Darcy's  flux 
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c)  Pore  water  velocity 
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t-t',  in  days 


 60%  water  uptake,  §^=1 

 60%  water  uptake,  §^="1 

 60%  constant  water  uptake 

—  —  No  water  uptake 


Figure  3.10:  Effects  of  the  water  uptake  and  of  the  cross-correlations  between  rainfall 
and  water  uptake  with  60  %  mean  water  uptake  scenario  at  the  initial  time  t'  =  60 
days  and  at  the  depth  z=z'=2  meters  on  the  temporal  behavior  of  the  autocorrelations 
of  a)  the  soil-water  content  b)  the  Darcy's  flux  and  c)  the  pore-water  velocity. 
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Figure  3.11:  Effects  of  the  water  uptake  and  of  the  cross-correlations  between  rainfall 
and  water  uptake  with  60  %  mean  water  uptake  scenario  at  the  initial  time  t'  =  60 
days  and  at  the  depth  z=z'=3  meters  on  the  temporal  behavior  of  the  autocorrelations 
of  a)  the  soil-water  content  b)  the  Darcy's  flux  and  c)  the  pore-water  velocity. 
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c)  Pore  water  velocity 


z-z',  in  meters 


 60%  water  uptake,  g^=l 

 60%  water  uptake,  §^=-l 

 60%  constant  water  uptake 

—  —  No  water  uptake 


Figure  3.12:  Effects  of  water  uptake  and  of  the  cross-correlations  between  rainfall 
and  water  uptake  with  60  %  mean  water  uptake  scenario  at  the  original  depth  z'=2 
meters  and  at  t=t'=  60  days  on  the  spatial  behavior  of  the  autocorrelations  of  a)  the 
soil-water  content  b)  the  Darcy's  flux  and  c)  the  pore-water  velocity. 
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a)  Soil  water  content  b)  Darcy's  flux 


z-z',  in  meters  z-z',  in  meters 


c)  Pore  water  velocity 


z-z',  in  meters 


 60%  water  uptake,  §^=1 

 60%  water  uptake,  g^=_l 

 60%  constant  water  uptake 

—  —  No  water  uptake 


Figure  3.13:  Effects  of  the  water  uptake  and  of  the  cross-correlations  between  rainfall 
and  water  uptake  with  60  %  mean  water  uptake  scenario  at  t=t'=  60  days  and  at  z'=3 
meters  on  the  spatial  behavior  of  the  autocorrelations  of  a)  the  soil-water  content  b) 
the  Darcy's  flux  and  c)  the  pore-water  velocity. 


CHAPTER  4 

UNSATURATED  TRANSIENT  STOCHASTIC  TRANSPORT  MODEL 

4.1  Introduction 

The  resident  concentration  and  the  solute  flux  concentration  are  the  two  pre- 
diction and  detection  modes  used  to  quantify  solute  mass  in  soil  and  groundwater 
(Kreft  and  Zuber,  1978).  The  first  expresses  the  mass  of  solute  per  unit  volume  of 
fluid  contained  in  an  elementary  volume  of  the  system  at  a  given  instant.  The  second 
expresses  the  mass  of  solute  per  unit  of  fluid  passing  through  a  cross-section  at  an 
elementary  time  interval.  The  choice  of  whether  to  use  one  or  the  other  is  governed 
by  their  underlying  assumptions.  The  prediction/detection  mode  of  resident  concen- 
tration measures  concentrations  in  a  particular  volume  assuming  that  the  positions 
of  the  solute  molecule  are  immobile  during  the  time  of  measurement  (Jury  and  Roth, 
1990),  whereas  in  the  second  prediction/detection  mode  the  molecules  are  measured 
as  they  move  past  some  boundary.  Although  it  is  difficult  in  practice  to  measure 
the  solute  flux  concentration  without  disturbing  the  flow  lines,  it  is  important  to 
realize  that  solute  flux  concentration  is  a  more  realistic  measurement  because  the  so- 
lute particles  are  mobile  in  underground  formations.  Therefore,  further  consideration 
should  be  given  to  this  latter  prediction/detection  mode  for  groundwater  monitoring 
program  and  field  scale  transport  studies. 

In  natural  geologic  formations  characterized  by  large  variability  of  hydraulic 
properties,  it  seems  that  the  solute  flux  concentration  is  the  more  suitable  predic- 
tion/detection mode.  Shapiro  and  Cvetkovic  (1988)  showed  that  the  random  struc- 
ture of  underground  formations  can  cause  a  significant  difference  between  the  two 
detection  modes.  Another  result  by  Dagan  et  al.  (1992)  also  supports  this  choice. 
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They  found  that  solute  concentration  pumped  out  of  a  well  penetrating  a  stratified 
formation  of  different  permeabilities  is  representative  of  the  solute  flux  concentration 
rather  than  the  resident  concentration. 

Recent  developments  of  the  solute  flux  theory  have  been  driven  by  works  in 
saturated  flow  conditions  (Simmons,  1982;  Shapiro  and  Cvetkovic,  1988;  Dagan  et 
al,  1992)  as  well  as  in  unsaturated  flow  conditions  (Dagan  and  Dresler,  1979;  Jury 
et  al,  1986;  Destouni,  1991,  1992a,  1992b,  1993;  Russo,  1993b,  Russo  et  al,  1994a, 
1994b)  and  in  soil-groundwater  systems  (Destouni  and  Graham,  1995).  Most  of 
the  aformentioned  works  simplified  the  flow  regime  and  treated  solute  flux  transport 
under  steady-state  conditions.  Russo  et  al.  (1994a,  1994b)  tested  the  reliability  of 
the  steady-state  Lagrangian-stochastic  analyses  of  vadose-zone  transport  developed 
by  Russo  (1993a,  1993b)  under  a  deterministic  transient  flow  regime  in  order  to 
predict  solute  flux  breakthrough  curve.  These  works  showed  that  the  steady-state 
Lagrangian-stochastic  analysis  fails  to  accurately  predict  the  solute  spreading  and 
breakthrough  under  transient  flow  regime  due  to  the  complicated  velocity  fields. 

Since  most  of  the  unsaturated  transport  models  assume  a  steady-state  flow 
regime,  they  can  not  be  applied  to  practical  situations  in  which  water  flow  in  the 
upper  part  of  the  soil  is  affected  by  transient  processes  such  as  irrigation,  water  up- 
take and  precipitation.  All  these  limitations  can  be  overcome  by  running  intensive 
Monte-Carlo  simulations  which  combine  statistical  techniques  to  generate  realiza- 
tions of  heterogeneous  underground  formations  with  efficient  numerical  schemes  that 
can  handle  the  complex  partial  differential  equations  of  unsaturated  flow.  However, 
besides  the  difficulty  of  obtaining  numerical  stability  and  convergence,  the  major 
drawback  of  this  type  of  approach  is  the  formidable  computational  effort  required  to 
solve  the  non-linear  unsaturated  flow  equation  in  large  and  heterogeneous  formations 
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under  transient  flow  regimes.  Therefore,  the  development  of  a  simplified  stochas- 
tic transport  model  for  heterogeneous  conditions  and  under  transient  flow  regime  is 
needed. 

The  general  objective  of  this  Chapter  is  to  develop  and  test  a  semi-analytical 
one-dimensional  stochastic  transport  model  for  inert  and  reactive  solutes  in  heteroge- 
neous formations  subject  to  a  transient  flow  regime.  The  general  method  predicts  the 
solute  flux  statistics  based  on  those  of  the  travel  time  using  a  Lagrangina  approach. 
First,  the  mean  and  the  variance  of  the  travel  time  distribution  are  determined  from 
the  statistics  of  the  pore-water  velocity  derived  in  the  previous  Chapters.  Then, 
the  travel  time  moments  are  used  in  an  assumed  travel  time  probability  distribution 
function  with  an  appropriate  impulse  response  function. 

4.2    Review  of  Solute  Flux  Theory 

The  one  dimensional  Eulerian  transport  equation  for  the  solute  flux  concen- 
tration in  the  unsaturated  zone  has  as  its  starting  point  the  well-known  mass  balance 
equation: 


d9(z,t)Cf(z,t)  +  dS(z,  t)  =  Q 
dt  dz 

where  S(z,t)  is  the  solute  flux  [^]  defined  by  q(z,t)Cf  (z,t)  in  terms  of  the 
solute  flux  concentration  Cf  [f]  or  by  q(z,t)  Cr{z,t)  -  0(z,t)D(z,t)^~^  in  terms 
of  the  resident  concentration  Cr  [^].  The  quantities  q  [£],  D  [4r],  and  6  denote 
the  Darcy's  flux,  the  longitudinal  local  dispersion  coefficient  and  the  volumetric  soil- 
water  content,  respectively.  Note  that  the  same  equation  is  also  valid  for  the  resident 
concentration  (Parker  and  Van  Genuchten,  1984)-  Assuming  that  q,  6  and  D  are 
constant,  the  solution  of  equation  (4.1)  with  the  following  initial  and  boundary  con- 
ditions, 
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Cf{z,0)  =  0 

Cf(z  =  0,t)  =  —S(t) 

Cf{z  =  +oo,0)  =  0 

is  given  by  Jury  and  Roth,  1990: 

,  M0z      _(*-Z')2  Mo 

CJ(z,t)  =  iDt    = — 1\v^) 

2qV7rDt3  Q 

where  M0  is  the  injected  mass  [M],  V  is  the  pore- water  velocity  [|?]  defined  by 
2.  Because  the  governing  partial  differential  equation  is  linear,  ^(V,  t)  [^}  is  equivalent 
to  the  impulse  response  function  of  the  linear  transport  system  to  an  instantaneous 
pulse  injection.  Therefore,  the  solute  flux  distribution  can  be  generalized  for  any 
given  temporal  distribution,  f(t),  of  total  mass  M0  introduced  into  the  system,  by 
the  convolution  integral: 

S{z,  t)  =  f  M0f(t  -  r)7(V,  r))dr  (4.2) 

Assuming  that  the  uncertainty  of  the  mass  balance  transport  equation  is  due 
to  the  velocity,  the  statistics  of  the  solute  flux  can  be  calculated  once  the  probability 
density  function  of  the  velocity  field  is  specified. 

The  same  formulation  is  obtained  for  the  solute  flux  using  the  Lagrangian 
framework  initially  developed  by  Dagan  et  al.  (1992)  and  then  generalized  by  Destouni 
and  Graham  (1995).  The  solute  flowing  through  the  porous  medium  is  conceived  as 
set  of  discrete  particles  of  indivisible  mass  moving  in  an  ensemble  of  streamtubes. 
The  flow  regime  is  assumed  to  be  steady  so  that  the  streamline  and  pathlines  coin- 
cide (Bear,  1972).  The  individual  solute  flux  within  one  streamtube  is  defined  by 
Destouni  and  Graham  (1995,  eq.  2)  as  a  convolution  integral, 
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S(z,  t,  a)  =  M0  f  f(t  -  r)7(T,  r,  a.)dr  (4.3) 
Jo 

where  M0  is  the  total  particle  mass  [M]  introduced  into  the  system,  /  is 
the  input  distribution  function  [^],  a  is  the  starting  position  of  the  particle,  7  is 
the  impulse  response  function  [£]  and  T  is  the  random  travel  time.  Contrasting 
equations  (4.2)  and  (4.3)  leads  to  the  conclusion  that  they  are  equivalent  in  the  sense 
that  the  random  velocity  of  the  former  corresponds  to  the  random  time  of  the  latter 
for  a  fixed  z.  Simmons  (1982)  has  discussed  and  developed  conversion  relationships 
between  these  two  approaches. 

Assuming  that  the  travel  time  is  random,  the  expected  value  and  the  variance 
of  the  solute  flux  are  given  by 

poo  pt 

(S(z,t))=  /     /  M0f(t-T)i{T,T)9l(T)dT  dT  (4.4) 
Jo  Jo 

a2(z,t)  =  (S2(z,t))  -  (S(z,t))2  (4.5) 

with 

poo  pt 

(S2(z,  t)}  =  Mq   /     /  [/(*-t)7(T,t)]2  9l(T)dr  dT 
Jo  Jo 

where  the  ()  stands  for  the  expected  value,  and  g\  is  the  one  particle  travel 
time  probability  density  function.  The  difficulty  in  this  Lagrangian  framework  is  the 
determination  of  an  appropriate  travel  time  distribution.  Experimental  measurements 
can  be  used  to  determine  the  travel  time  distribution,  but  with  this  approach  the 
results  are  site-specific  and  can  not  be  extrapolated.  Therefore,  the  alternative  is  to 
derive  the  travel  time  ensemble  moments  as  a  function  of  the  formation  properties 
and  weather  statistics  that  are  more  easily  quantified. 
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4.3    Analytical  Development 
4.3.1    Travel  Time  Statistics  for  Non-reactive  Solute 

The  dimensionless  travel  time,  T,  of  a  particle  moving  through  a  vertical 
porous  medium  from  depth  z=0  to  z=Z  is  described  by, 


dz 


r  t 

T(Z)  =  /  -j- 
Jo  v(z, 


;*(*)) 

where  v(z,t(z))  represents  the  local  dimensionless  pore-water  velocity  defined 
as  v(z,t(z))  =  fffiz))  wnere  Q(z^(z))  and  Q(zit(z))  are  tne  local  dimensionless 
Darcy's  flux  and  the  volumetric  soil- water  content,  respectively.  Using  a  Taylor's 
series  expansion  around  the  zero  order  pore-water  velocity  i>0oo  (Mood  et  al.,  1974), 
the  travel  time  becomes 

1       (v(z,t(z)) -v000) 


-f 

Jo 


Vqoo  v$00 


+ 


dz 


where  i>ooo  denotes  the  zero  order  mean  pore-water  velocity  defined  previously 
in  Chapters  2  and  3  and  all  scaling  factors  required  to  non-dimensionalize  these 
equations  are  as  defined  in  Chapters  2  and  3.  Evaluating  the  mean  travel  time  at  the 
zero  and  first  order  leads  to: 

CO  -  i  (4-6) 

The  travel  time  variance  of  the  inert  solute  is  approximated  by  the  double 
spatial  integral  of  the  second  order  pore-water  velocity  autocovariance  determined  in 
equations  (2.51)  and  (3.31)  (Dagan  et  al,  1992,  eq.  14): 

4iZ)={[Z±-  ("(»■«(»)) fZJ__  Wj^Wh^oo)^) 
Jo    vooo  Vqoo  Jo    vooo  v000 

=~i-  /     /    Pw(z,t(z);z',t(z'))dzdz'  (4.7) 
%)0  Jo  Jo 
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The  time  t  corresponding  to  a  given  position  z  of  the  particle  is  random  and 
therefore  unknown.  Thus,  it  is  assumed  that  this  time  can  be  approximated  by  its 
first  order  mean,  t  ~  (T)  ~  —  rather  that  at  the  random  instantaneous  velocity 
(Dagan,  1984;  Zhang  and  Neuman,  1996).  This  leads  to  the  following  second  order 
travel  time  variance: 

o2T{Z)~  I*  [  Pvv(z,—;z',—)dzdz'  (4.8) 
uooo  Jo    Jo  vooo  ^000 

Figure  4.1  illustrates  the  evolution  of  the  integrand  of  equation  (4.8)  in  a 
heterogeneous  fine  sandy  soil  subject  to  a  rainfall  pattern  representative  of  Tampa, 
Florida,  when  the  times  of  observation  are  adjusted  to  the  depth  of  observation 
through  the  relationship  t  ~  Note  that  the  integrand  is  not  stationary  and 

becomes  more  dispersed  when  depth  increases. 

The  dimensional  mean  travel  time  is  obtained  in  multiplying  equation  (4.6) 
by       for  heterogeneous  soil  and  ^  for  homogeneous  soils.  The  dimensional  travel 

G  3 

time  variance  is  obtained  by  multiplying  equation  (4.8)  by  (^-)2  for  heterogeneous 

G 

soil  and  (pi)2  for  homogeneous  soils. 

4.3.2    Travel  Time  Statistics  for  Reactive  solute 

The  dimensionless  travel  time,  Tr,  of  a  reactive  solute  moving  through  a  ver- 
tical porous  medium  from  depth  z=0  to  z=Z  is  described  by, 


T, 


where  vr  is  the  instantaneous  dimensionless  pore-water  velocity  of  the  reactive 
solute  defined  as: 


Vr(z,t)  =  Q,  .  ■  (4.10) 

9{z,t)  +pKd(z)  v  ' 
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with  p  and  Kd  are  the  uniform  bulk  density  of  the  soil  and  the  spatially  variable 
linear  equilibrium  sorption  coefficient,  respectively.  Note  that  the  product  of  pKd(z) 
is  dimensionless.  Substituting  equation  (4.10)  into  equation  (4.9)  after  expanding  Kd 
into  the  sum  of  mean  and  a  perturbation,  Kd(z)  =  (Kd)  +  uKdl(z)  leads  to: 

^J^JM,  (4,1) 
J0    v(z,t{z))     J0  q(z,t(z)) 

where  Kdl  denotes  a  dimensionless  random  space  function  with  a  zero  mean 
and  an  unit  standard  deviation,  and  u  =  OKd  represents  the  standard  deviation  of 
Kd-  Note  that  the  reactive  travel  time  is  made  of  the  sum  of  an  inert  component 
corresponding  to  the  travel  time  developed  in  the  previous  section  and  of  a  reactive 
component  corresponding  to  the  effect  of  the  spatially  random  linear  equilibrium 
adsorption  coefficient  Kd. 

Evaluating  the  dimensionless  mean  travel  time  of  reactive  solute  at  the  zero 
and  first  order  leads  to: 

(Tr(Z))  =  —  +  ^  (4.12) 
^ooo  9ooo 

Evaluating  the  approximate  travel  time  variance  requires  substituting  equation 
(3.27)  and  (3.30)  into  equation  (4.11)  for  the  travel  time  at  location  z  and  time  t, 
multiplying  by  a  different  travel  time  at  location  z'  and  time  t',  taking  the  expected 
value  of  the  resulting  equation,  and  retaining  terms  up  to  the  second  order: 
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^000 
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(4.13) 


where  urooo  =  e00o+p{Kd)  ls  ^e  ^rst  orc^er  pore-water  velocity  of  the  reactive 
solute.  Note  that  the  pore- water  velocity  covariances  in  equation  (4.13)  are  evaluated 
by  adjusting  the  times  and  depth  of  observation  around  the  first  order  mean  pore- 
water  velocity  of  the  reactive  solute  (Bellin  et  al,  1993).  Note  that  the  travel  time 
variance  is  also  composed  of  a  inert  component  corresponding  to  that  derived  in 
the  previous  section  and  a  reactive  component  that  combines  the  effect  of  rainfall, 
water  uptake,  and  linear  equilibrium  adsorption.  PeomKdl  (z,  z')  is  the  cross-covariance 
between  the  spatial  component  of  the  soil-water  content  and  the  perturbation  of 
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the  linear  equilibrium  adsorption  coefficient.  Taking  the  expected  value  of  equation 
(2.17)  (which  is  equivalent  to  the  governing  equation  of  0Oio)  and  its  associated  initial 
and  boundary  conditions  after  postmultiplying  by  Kdl  (z) ,  leads  to  the  steady-state 
governing  equation  of  Pe010Kd  {z,  z')  (valid  away  from  the  initial  condition): 


V(6000) 


d2Pe010Kdl  (z1;  z)  dfC 


dz 


1 2 
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dPe010Kdl  {z'\  z)  _  d_ 
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=  —  JC(900o)PVKdl  (z\  z) 
z>=0;z  (4.15) 


=  — /C(^ooo)-PrjA'<il  (z',  z) 


z'=+<x;z 


(4.16) 


where  PVKdl  {z1,  z)  denotes  the  cross-covariance  between  the  dimensionless  per- 
turbations of  the  saturated  hydraulic  conductivity  and  linear  equilibrium  adsorption 
coefficient.  The  forcing  terms  in  equations  (4.14)-(4.16)  are  identical.  Thus,  if  PnKdi 
is  assumed  stationary  over  space,  an  infinite  domain  solution  can  be  found  that  sat- 
isfies these  equations.  Taking  the  Fourier  transform  of  equation  (4.14)  leads  to  the 
cross-spectrum  between  0Oio  and  , 


Se0ioKdl  (h)  — 


i  £(0ooo)  SVKdl(h) 
-^(^ooo)^i  -  i 


dK\ 
dO  1 0ooo 


(4.17) 


The  cross-covariance,  Pe010Kd  is  found  by  taking  the  inverse  Fourier  transform 
of  Sg010Kdl  after  SVKdl  has  been  specified.  In  order  to  calculate  the  travel  time  variance 
of  inert  or  reactive  solute,  the  temporal  components,  Pe1009100,  P0loooOoi>  and  Pe0oie0oi 
and  the  spatial  component  P9010e010  of  the  soil-water  content  autocovariance  were 
determined  following  the  procedure  described  in  Chapters  2  and  3.  This  involves 
using  a  Fast  Fourier  Transform  code  for  the  inversion  of  the  double  Laplace  transform 
of  each  temporal  component  and  using  the  inverse  Fourier  Transform  for  the  spatial 
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component.  Then,  the  travel  time  variance  was  calculated  numerically  using  the 
trapezoidal  rule. 

The  dimensionless  travel  time  mean  and  variance  of  the  reactive  solute  in  a 
heterogeneous  soil  are  dimensionalized  by  dividing  both  the  inert  and  reactive  com- 
ponents of  equations  (4.12)  and  (4.13)  by  K*G  and  JC*G2,  respectively.  Similarly,  the 
dimensionless  travel  time  mean  and  variance  of  the  reactive  solute  in  a  homogeneous 
soil  with  respect  to  saturated  hydraulic  conductivity  are  dimensionalized  by  i)  setting 
/  =  0,  and  ii)  dividing  the  inert  component  of  equations  (4.12)  and  (4.13)  by  /C*  and 
JC*2,  respectively.  Consideration  of  the  water  uptake  implies  replacing  0Ooo  by  0Ooo- 

4.3.3    Solute  Flux  Statistics 

A  stochastic  advective-dispersive  streamtube  transport  model  predicts  the 
mean  and  variance  of  solute  flux  by  using  equations  (4.4)  and  (4.5)  that  require 
specifying  a  impulse  response  function  and  a  travel  time  probability  distribution. 
For  illustrative  purposes,  a  lognormal  advective-dispersive  impulse  response  function, 
j(T,  t),  was  used  for  a  conservative  solute  applied  as  a  instantaneous  pulse  at  the  soil 
surface: 


1  (In(t)-ln(T))2 

7(T,t)  =  -=— c  (4.18) 

Over  depth,  the  response  impulse  function  is  characterized  by  the  following 
moments  taken  from  Goltz  and  Roberts  (1987): 


2 

9ooo  <7ooo 
2         {ZOooo)2  .  6DZ9lQ0     12D  8q00  2 

°L     =   2  '  3  '  4  

%00  %00  %00 

where  T  is  the  local  arrival  time  of  the  center  of  mass  determined  by  the  local 
zero  order  mean  soil- water  content  and  its  corresponding  mean  Darcy  flux,  g0oo,  D 
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is  the  local  deterministic  hydrodynamic  dispersion  coefficient  within  the  streamtube, 
o\  is  the  local  spread  around  T  caused  by  D,  and  o\  =  ln((^)2  +  1)  is  the  spatial 
variance  associated  with  the  impulse  response  function. 

The  lognormal  distribution  has  been  hypothesized  for  the  travel  time  proba- 
bility density  function.  Early  investigations  have  substantiated  this  choice  (Jury  and 
Roth,  1990  and  Russo  et  al,  1994a).  Therefore,  in  this  work  the  single  particle  travel 
time  distribution,  gi{T),  is  represented  by, 


(in(T)-(ln(T)))2 
7i 


9l(T)  =  e       HTt  (4.19) 

V 1-KOlnTl 


where  afnT  =  Ln{-^  +  1)  and  (ln{T))  =  ln{(T))  -  Assuming  a  dirac 
input  function  at  the  soil  surface  and  substituting  equations  (4.18)  and  (4.19)  into 
equation  (4.4)  leads  to  the  mean  solute  flux  passing  a  given  depth  z, 


rOO  I  (In(t)-ln(T))2  1  _  (In(T)-(Jn(T)))2 

(S(z,t))  =  M0        -=^-e  e  dT 

Jo    V2TraDt  V ZnainTi 

After  making  the  changes  of  variables  x  =  Int  and  y  =  InT,  the  analytical 

solution  of  the  previous  equation  is  found  in  Abramovitz  and  Stegun  (1970): 


<SM>  =  e*^  (4.20) 


where  a,  =  h  =  -Wnp^mt)  and     =  rtw'+^fyint+of^nt)  ^ 

Using  the  aformentioned  changes  of  variable,  evaluating  the  solute  flux  vari- 
ance requires  first  evaluating  (S(z,  t)2)  , 
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where  a\  =  <&J§^,  b[  =  -Wff^nr),  and 

,  _  („%  (lnT)2+Aln  t<T2n^T+2(fn  t)2<rfnT) 

Therefore,  the  solute  flux  variance  within  one  streamtube  is  established  as: 

a2s(t)   =   (S(z,t)2)  -  (S(z,t))2 

poo 

=  MS  J  MT,t))2gi(T)dT-(s(t)f 

(27r)faf,<TjnT  V  ai  2naDolnT  \j  ax 

where  ai,  61,01,^,61,  and c'x  are  as  defined  above. 

4.4    Illustration  of  Theoretical  Results 

The  effect  of  saturated  hydraulic  conductivity,  water  uptake  and  the  equilib- 
rium sorption  coefficient  on  the  statistics  of  the  travel  time  and  solute  flux  will  be 
investigated  for  inert  and  reactive  solutes  moving  in  a  fine  sandy  soil  subject  to  ran- 
dom rainfall  pattern  representative  of  Tampa,  Florida.  The  local  relative  hydraulic 
conductivity,  K{0),  and  the  soil-moisture  diffusivity  V  are  described  by  the  Brooks 
and  Corey  model  (Brooks  and  Corey,  1964)-  As  an  illustration,  predictions  of  reactive 
solute  transport  are  made  for  the  pesticide  Napropamide  (Jury  and  Gruber,  1989). 
The  effect  of  water  uptake  is  undertaken  for  the  60%  mean  uptake  scenario  presented 
in  the  previous  Chapter.  The  statistics  of  the  random  parameters  are  summarized  in 
Table  4.1  while  the  deterministic  properties  of  the  soil  are  presented  in  Table  4.2. 

Inert  Solute 

Travel  time  and  solute  flux  statistics  calculated  for  random  rainfall  and  uniform  sat- 
urated hydraulic  conductivity  and  random  rainfall  and  heterogeneous  saturated  hy- 
draulic conductivity  are  presented  on  Figures  4.2a,  b,  and  c.  It  is  obvious  from  these 
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figures  that  the  saturated  hydraulic  conductivity  variability  does  not  have  a  signifi- 
cant impact  on  the  one-dimensional  transport  at  the  selected  depths  for  a  variance 
of  lnK*s  =  1. 

Figures  4.3  and  4.4  demonstrate  the  effect  of  water  uptake  at  two  different 
depths  and  for  two  different  levels  of  cross-correlation,  psx,  between  rainfall  and  water 
uptake,  on  the  transport  statistics.  The  effect  of  the  water  uptake  has  been  evaluated 
with  the  60%  uptake  scenario  presented  in  the  previous  Chapter.  Compared  to  no 
water  sink,  one  can  observe  that  water  uptake  with  any  level  of  psx  causes  i)  a  delay 
in  the  zero  order  mean  travel  time  to  a  given  depth,  and  ii)  an  increase  of  the  second 
order  variance  of  the  travel  time.  As  a  result,  the  mean  and  the  standard  deviation  of 
the  solute  flux  breakthrough  curves  will  exhibit  lower  peaks  and  larger  spreading.  This 
phenomenon  is  enhanced  as  the  cross-correlation  between  water  uptake  and  rainfall 
becomes  negative.  From  the  effect  of  water  uptake  on  the  flow  statistics  discussed  in 
Chapter  3,  such  results  are  expected.  As  noted  earlier,  the  effect  of  the  water  uptake 
is  i)  to  reduce  the  mean  pore- water  velocity,  ii)  to  slightly  decrease  the  variance  of  the 
pore-water  velocity,  and  iii)  to  increase  significantly  the  temporal  correlation  scale  of 
the  pore-water  velocity.  The  delay  of  the  mean  travel  time  is  due  to  the  reduction  of 
the  mean  pore-water  velocity.  The  increase  of  the  travel  time  variance  is  due  to  the 
double  integration  of  a  pore-water  velocity  autocovariance  which  is  correlated  over 
larger  times. 

These  theoretical  results  do  not  agree  with  the  findings  of  Russo  et  al.  (1998). 
Performing  a  single,  three-dimensional,  numerical  simulation  in  a  heterogeneous,  un- 
saturated soil  under  transient  flow  regime,  Russo  et  al.  examined  the  water  uptake 
effect  on  solute  flux  concentration  and  on  the  travel  time  probability  density  func- 
tion. Using  a  distribution  of  normalized  solute  flux  concentration,  Russo  et  al.  (1998) 
show  that  the  effects  of  water  uptake  are  to  reduce  the  skewing  and  to  increase  the 
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peak  of  the  travel  time  distribution.  However,  they  observed  that  the  standard  de- 
viation of  the  normalized  solute  flux  concentration  (their  Figure  16)  increases  with 
the  presence  of  water  uptake.  In  the  Lagrangian  approach  to  transport,  these  results 
appear  to  be  contradictory.  On  one  hand,  the  water  uptake  reduces  the  skewing  of 
the  travel  time  distribution  which  implies  a  decrease  of  the  travel  time  variance.  On 
the  other  hand,  the  effect  of  water  uptake  is  to  increase  the  standard  deviation  of 
the  normalized  solute  flux  concentration,  which  implies  an  increase  of  the  travel  time 
variance.  Several  modeling  differences  exist  between  the  previous  mentioned  work 
and  the  current  model.  Their  single  numerical  simulation  may  not  be  extrapolated  to 
ensemble  statistics.  Their  water  uptake  flux  was  a  function  of  the  soil-water  content 
whereas  in  the  current  stochastic  transport  model  the  water  uptake  is  a  function  of 
an  exponential  decaying  function  multiplied  by  a  temporally  random  water  sink.  The 
results  of  Russo  suggest  that  it  would  be  of  interest  to  run  multiple  simulations  with 
random  rainfall  and  water  uptake,  and  to  examine  the  temporal  as  well  as  the  spatial 
correlation  scales  of  the  pore-water  velocity  under  such  conditions. 

Reactive  Solute 

The  effect  of  equilibrium  sorption  coefficient  on  transport  statistics  were  simulated 
in  a  heterogeneous  fine  sandy  soil  subject  to  the  Tampa  random  rainfall.  The  cross- 
covariance  between  the  linear  equilibrium  sorption  coefficient  and  the  saturated  hy- 
draulic conductivity  was  derived  from  assuming  that  the  linear  sorption  coefficient 
is  a  weighted  combination  of  two  fields,  one  correlated  to  the  saturated  hydraulic 
conductivity  and  one  uncorrected  (Gelhar,  1993),  that  is 

Kdl  =  PvKdl  T)  +  ^l-p2nKdi  V 

where  v  is  a  random  zero-mean  spatially  correlated  process  that  is  uncorrelated 
with  77,  and  pnKdl  is  a  weighting  parameter  (—1  <  pVKdl  <  1)  that  determines  the 
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correlation  between  r]  and  K^.  The  covariances  used  for  the  reactive  solute  are 
presented  in  Table  4.3. 

Inverting  equation  (4.17)  with  Mathematica  (Wolfram,  1991)  after  substitu- 
tion  of  Sr,Kdl  =  ^(1+A2fc5)   leads  to: 

(P(gooo)e^  -  Vf  L000  e~  ) 

where  £  =  \z  —  z'\.  The  form  of  the  lognormal  impulse  response  function  for 
the  reactive  solute  were  assumed  to  be  identical  to  these  of  the  inert  solute.  The 
effect  of  equilibrium  sorption  on  transport  statistics  is  shown  on  Figures  4.5  and 
4.6  in  a  heterogeneous  fine  sandy  soil  for  two  different  cross-correlations  between 
saturated  hydraulic  conductivity  and  equilibrium  sorption  coefficient,  Pr)Kdl  =  1  and 
pVKd  =  —  1.  The  effect  of  the  equilibrium  sorption  on  the  travel  time  are  i)  to  increase 
the  mean  and  also  the  variance  of  the  travel  time  compared  to  the  inert  solute,  ii) 
to  delay  and  reduce  the  peak  of  the  breakthrough  curves  mean  and  variance  of  the 
solute  flux,  and  iii)  to  enhance  spreading  of  the  solute  flux  statistics.  Note  that  the 
sign  of  pVKdl  has  negligible  effect  on  the  statistics  of  travel  time  and  solute  flux  under 
a  transient  flow  regime. 

4.5    Validation  Against  Monte-Carlo  Simulations 

In  the  absence  of  a  large  and  reliable  collection  of  field  data,  a  set  of  intensive 
Monte-Carlo  experiments  was  implemented  in  order  to  test  the  performance  of  the 
stochastic  advective-dispersive  streamtube  transport  model. 

4.5.1    Monte  Carlo  Simulations 

The  domain  of  simulation  represents  a  vertical  cross  section  of  of  a  30  m  wide 
by  10  m  deep  soil  overlying  a  horizontal  water  table  (Fig.  4.7).  Rainfall  enters  the 
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system  uniformly  over  the  entire  land  surface.  No  horizontal  flow  occurs  on  either 
side  of  the  domain.  A  modified  version  of  the  variably  saturated  flow  and  transport 
code,  VS2DT  (Lappala  et  al,  1993,  Healy,  1990),  was  used  to  simulate  the  infiltra- 
tion experiment.  This  code  solves  the  head-based  Richards  equation  for  flow  and 
the  advection-dispersion  equation  for  transport.  The  Brooks  and  Corey  model  was 
selected  to  model  the  soil  moisture  characteristic  curves.  The  code  was  modified 
in  order  to  accept  spatially  random  inputs  of  saturated  hydraulic  conductivity  and 
temporally  random  time  series  of  rainfall  fluxes.  Appropriate  spatial  and  temporal 
discretization  was  adopted  to  insure  convergence,  stability  and  accuracy  of  the  nu- 
merical scheme.  Vertical  and  horizontal  grid  sizes  of  .25  m  and  a  maximum  time 
step  of  .01  days  satisfy  the  Peclet,  Neuman  and  the  Courant  number  conditions.  The 
initial  condition  soil-water  content  was  established  using  the  Brooks  and  Corey  equi- 
librium pressure  profile  above  the  horizontal  water  table.  Rainfall  was  treated  as  a 
temporally  random  boundary  condition  at  the  land  surface. 

Intensive  Monte  Carlo  simulations  including  150  runs  each  were  performed  for 
two  rainfall  patterns  representative  of  Tampa  (Florida)  and  Olympia  (Washington 
State)  and  two  soil  types  for  which  the  saturated  hydraulic  conductivity  was  assumed 
to  be  either  uniform  or  spatially  variable.  The  random  saturated  conductivity  were 
generated  using  the  turning  band  algorithm,  originally  developed  by  Journel  and 
Huijbregts,  (1978)  and  then  modified  by  Akpoji  (1993).  From  Gelhar  (1993,  p.  291), 
an  averaged  value  of  1  was  set  for  the  variance  of  InfC*,  afnlc,  =  /  =  1.  The  rainfall 
series  were  generated  using  WGEN,  a  weather  generator  developed  by  Richardson 
and  Wright  (1984).  This  computer  simulation  code  generates  daily  rainfall  series 
which  reproduces  the  statistics  of  the  rainfall  from  a  specified  region.  WGEN  uses  a 
first-order  Markov  process  to  represent  the  relative  occurences  of  wet  day  (P\v)  and 
dry  days  (PD)  and  a  gamma  distribution  to  represent  the  daily  volume  of  the  rainfall. 
Statistical  parameters  of  the  rainfall  series  and  saturated  hydraulic  conductivity  fields 
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are  summarized  in  Table  4.4.  The  simulations  were  run  for  a  fine  sandy  soil  and 
a  sandy  loam  for  which  the  deterministic  properties  of  the  moisture  release  curves 
are  presented  in  Table  4.5.  In  order  to  compare  the  analytical  development  with 
the  numerical  experiments,  a  negative  exponential  function  was  fit  to  the  sample 
covariance  function  produced  by  WGEN  for  each  weather  scenario.  Note  that  annual 
average  rainfall  statistics  were  implemented  to  maintain  stationarity  of  the  rainfall 
time  series.  Values  of  solute  flux  were  recorded  at  one  node  located  in  the  middle  of 
four  control  planes  CP1,  CP2,  CP3,  and  CP4  and  relevant  statistics  calculated  for 
each  of  the  400  days  and  over  all  the  runs.  For  the  fine  sandy  soil,  the  four  control 
planes  are  located  at  .6125,  2.125,  4.125,  and  6.125  m  from  the  land  surface,  whereas 
for  the  sandy  loam  soil  the  four  controls  planes  are  located  at  .6125,  1.125,  2.125,  and 
3.125  m. 

4.5.2    Results  and  Discussion 

Figures  4.8-4.13  show  how  the  theoretical  (or  predicted)  mean  and  coefficient 
of  variation  of  the  solute  flux  produced  by  the  stochastic  advective-dispersive  un- 
saturated transport  model  compare  to  their  numerical  (or  measured)  counterparts. 
Figures  4.8a,  4.10a,  and  4.12a  demonstrate  the  good  performance,  except  at  shallow 
depths,  of  the  model  in  predicting  the  mean  solute  flux  in  fine  and  coarse  textured 
soils  subject  to  two  different  rainfall  patterns.  However,  at  shallow  depths,  the  model 
seems  to  systematically  overpredict  the  peak  mean  solute  flux.  This  discrepancy  may 
be  due  to  the  fact  that  at  shallow  depths  the  assumption  of  a  small  variability  of  the 
soil-water  content  required  by  the  stochastic  unsaturated  flow  model  is  violated.  In 
fact  at  shallow  depths,  the  rainfall  variability  causes  the  more  non-linearity  of  the 
flow  equation  solved  by  VS2DT  than  at  deeper  depths.  However,  it  is  worth  not- 
ing that  the  quality  of  the  predictions  improves  over  depth.  Compared  to  previous 
attempts  of  predicting  the  mean  solute  flux  as  illustrated  by  the  work  of  Russo  et 
al. (1994b),  the  current  model  demonstrates  a  significant  improvement  in  predicting 


148 


mean  solute  flux  under  transient  flow  regime.  Figures  4.8b-4.13b  show  that  the  min- 
imum values  of  the  standard  deviation  are  systematically  underpredicted  as  a  result 
of  first  order  approximation  but  that  the  stochastic  advective-dispersive  unsaturated 
transport  model  reproduces  the  general  evolution  of  the  coefficient  of  variation.  The 
model  describes  reasonably  well  the  concave  shape  that  reaches  a  minimum  value  at 
the  temporal  center  of  mass  of  the  mean  solute  flux  breakthrough  curve. 

Figures  4.8b,  4.9b,  4.10b,  and  4.11b  show  that  the  standard  deviations  are 
better  predicted  in  uniform  soils  than  heterogeneous  soils  and  in  coarse  textured  soils 
than  fine  textured  soils.  Soil  heterogeneity  creates  a  two-dimensional  flow  that  is  not 
modeled  by  the  one-dimensional  stochastic  unsaturated  flow  model,  as  mentioned 
in  Chapter  2.  Furthermore,  the  high  storage  capacity  of  a  fine  textured  soil  like 
the  sandy  loam  soil  tends  to  exhibit  more  a  two  dimensional  flow  compared  to  a 
coarse  textured  soil  like  the  fine  sandy  soil.  Therefore,  the  one-dimensional  stochastic 
advective-dispersive  unsaturated  transport  model  shows  its  limitations  in  predicting 
two-dimensional  flow  due  to  soil  heterogeneity  in  fine  textured  soils.  In  Figures  4.10b, 
4.11b,  4.12b,  and  4.13b,  the  quality  of  the  standard  deviation  predictions  is  shown 
to  be  better  in  the  fine  sandy  soil  under  the  Olympia  rainfall  pattern  than  in  the 
same  soil  under  the  Tampa  rainfall  pattern.  The  effect  of  the  rainfall  variability 
on  the  transport  can  be  explained  by  recalling  the  effect  of  rainfall  variability  on  the 
flow  statistics.  These  limitations  of  the  stochastic  unsaturated  flow  have  already  been 
discussed  in  Chapter  2  and  are  only  summarized  here.  The  linearized  functions  of  soil 
diffusivity  and  hydraulic  conductivity  about  the  zero  order  mean  soil-water  content 
do  not  adequately  model  the  non-linear  behavior  of  these  latter  functions,  especially 
at  shallow  depths  where  the  soil-water  content  varies  considerably  over  time  as  a 
result  of  high  variance  of  rainfall  and  long  time  intervals  between  rainfall  events. 

The  discrepancies  of  the  transport  statistics  between  the  numerical  experi- 
ments and  the  semi-analytical  stochastic  transport  have  been  explained  so  far  as 
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a  consequence  of  the  water  flow  by  looking  at  the  assumption  underlying  the  one- 
dimensional  flow.  The  assumptions  underlying  the  transfer  function  approach  can 
also  shed  some  light  on  such  results.  The  solute  flux  statistics  are  calculated  assum- 
ing a  linear  time-invariant  impulse  response  function  for  the  transport  process.  The 
reliability  of  such  model  lies  in  the  careful  choice  of  the  impulse  response  function. 
Once  that  function  has  been  chosen,  it  is  not  supposed  to  change  over  time.  Con- 
trasting this  assumption  with  the  Monte-Carlo  simulation  points  out  the  weakness 
of  assuming  a  constant  impulse  response  function  under  a  transient  flow  regime.  As 
a  consequence  of  the  flow  conditions  changing  over  time,  the  impulse  response  func- 
tion should  also  be  adjusted  over  the  time  so  that  it  accurately  reflects  the  response 
of  the  system  at  any  time.  In  fact,  assuming  a  constant  impulse  response  function 
is  equivalent  to  assuming  a  "mean"  impulse  response  function  that  is  not  likely  to 
adequately  apply  to  any  realization  of  the  rainfall  series  or  soil  profile.  A  similar 
problem  was  noted  by  Jury  and  Roth  (1990,  p.  108)  in  their  steady-state  stochastic 
analysis  of  solute  flux.  Jury  and  Roth  noted  that  when  the  net  applied  water  changes 
over  time,  the  amount  of  water  required  to  push  the  solute  through  the  soil  profile  is 
controlled  by  the  initial  storage.  Therefore,  they  attempted  to  overcome  this  problem 
by  developing  a  transient  solute  flux  probability  density  function  that  accounts  for 
the  initial  water  storage  of  the  soil  at  the  time  that  the  solute  is  applied. 

4.6  Summary 

This  Chapter  describes  the  development  and  testing  of  a  one-dimensional 
stochastic  advective-dispersive  unsaturated  transport  model  under  transient  flow  con- 
ditions. Following  a  Lagrangian  approach,  the  random  solute  flux  is  modeled  as  the 
convolution  of  a  time-invariant  lognormal  advective-dispersive  impulse  response  func- 
tion and  a  Dirac  input  of  solute  at  the  soil  surface.  The  statistics  of  the  solute  flux 
are  determined  assuming  a  lognormal  probability  density  function  for  the  travel  time 
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for  which  the  first  and  second  moments  are  estimated  from  the  ensemble  moments 
of  the  rainfall  series,  the  water  uptake,  the  equilibrium  sorption  coefficient,  and  the 
saturated  hydraulic  conductivity,  using  the  one-dimensional  Richards  equation.  This 
stochastic  transport  model  was  then  compared  to  intensive  Monte-Carlo  simulations 
for  two  rainfall  patterns  representative  of  Florida  and  of  Washington  State  and  two 
soil  types  of  contrasting  textures. 

The  main  findings  of  this  Chapter  are  summarized  below: 

-  Theoretical  as  well  as  numerical  experiments  have  demonstrated  that  for 
shallow  soils  in  humid  climates  the  spatial  variability  of  the  saturated  hydraulic  con- 
ductivity has  a  much  smaller  effect  on  vadose  zone  solute  transport  predictions  than 
the  temporal  variability  of  rainfall. 

-  Theoretical  developments  suggest  that  the  effect  of  water  uptake  is  to  increase 
the  mean  and  the  variance  of  the  travel  time.  Thereby,  water  uptake  reduces  the  peak 
and  enhances  spreading  of  the  breakthrough  curve  of  the  mean  solute  flux  as  well  as 
that  of  the  standard  deviation. 

-  The  incorporation  of  a  spatially  random  retardation  factor  into  the  transport 
formulation  increases  the  mean  and  the  variance  of  the  solute  travel  time  and  lowers 
the  peak  and  enhances  the  spreading  of  the  breakthrough  curves  of  the  mean  and 
standard  deviation  of  the  solute  flux. 

-  Except  at  shallow  depths,  the  stochastic  transport  model  accurately  predicts 
the  mean  solute  flux  breakthrough  curves  for  all  rainfall  patterns  and  soil  types  ex- 
amined by  the  Monte-Carlo  simulations.  Near  the  land  surface,  the  accuracy  of  the 
linearization  of  the  soil  hydraulic  properties  around  the  zero  order  mean  soil-water 
content  deteriorates  because  the  rainfall  pattern  causes  the  soil- water  content  to  vary 
significantly  over  time. 

-  The  stochastic  dispersive-advective  transport  model  is  able  to  reproduce  the 
general  concave  shape  and  magnitude  of  the  solute  flux  uncertainty  for  the  soil  and 
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climates  combinations  studied  here.  However,  the  model  does  systematically  under- 
predict  the  smallest  values  of  the  coefficient  of  variation.  The  solute  flux  coefficients  of 
variation  are  better  predicted  under  uniform  soil  conditions  than  under  heterogeneous 
conditions. 

-  The  Fast  Fourier  Transform  provides  good  approximations  of  the  solute  flux 
statistics  at  significant  computational  savings.  It  can  be  roughly  estimated  that 
computing  the  solute  flux  statistics  for  each  Monte-Carlo  run  generally  needs  N2  order 
operations  per  time  step  for  N  discrete  nodes  for  a  two-dimensional  problem.  The 
numerical  Fast  Fourier  Transform  involves  only  on  the  order  oiA*Log2N  operations  for 
N  nodes.  The  CPU  time  used  for  these  two-dimensional  randomly  transient  boundary 
problems  using  a  Sun  Ultra  workstation  is  about  three  hours  for  the  Fast  Fourier 
Transform  approximation  versus  twenty-five  days  using  a  set  of  150  Monte-Carlo 
runs. 
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Table  4.1:  Statistical  properties  of  the  rainfall  series,  the  log-saturated  hydraulic 
conductivity,  the  water  uptake,  and  linear  equilibrium  adsorption  coefficient. 


Variable 

lnK*t 

e* 

Mean 
Variance 

0.0038  m/d 
9.81  x  lO-W/d2 

1.61  m/d 
1 

.0023  m/d 
4  x  10-*m2/cP 

.08  dm3 /kg 
.0064  dm6 /kg2 

Correlation 

length 
Coef.  of 

0.6  d 

1  m 

.6  d 

1  m 

variation 

2.53 

.62 

.87 

1 

Table  4.2:  Deterministic  hydraulic  properties  of  the  fine  sandy  soil. 


Property 

Value 

0s 

0.35 

9r 

0.063 

K 

-.82  m 

X 

3.7 

p* 

1.2  kg/ 'dm3 

Table  4.3:  Covariances  used  for  the  statistics  of  the  reactive  solute  and  their  Fourier 
transforms. 


Covariance 

Expressions 

Fourier  transform 

p 
p 

P lKd1 
PKdlKdl 

e  x» 

_i£i 

e  Xr> 

PvKdl  °»m 

°w  -  ir(l+AJ*f) 

SKdl  Kdl  =  f%Kdl        +  Pr,/^ 
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Table  4.4:  Ensemble  statistics  of  the  two  rainfall  scenari  combined  with  the  two 
soils,  SI:  Tampa  scenario  x  fine  sandy  soil,  S2:  Olympia  scenario  x  fine  sandy 
soil,  S3:  Olympia  scenario  x  sandy  loam  soil.  P^W/D]  and  P^VF/W]  stand  for  the 
probabilities  of  a  wet  day  on  day  i  given  a  dry  day  on  day  i-1,  and  the  probability  of 
a  wet  day  on  day  i  given  a  wet  day  on  day  i-1,  respectively. 


SI 

S2 

S3 

Pi[W/W] 

.453 

.810 

.816 

Pi[W/D] 

.240 

.452 

.452 

Pw 

.30 

.71 

.71 

(Qo)  (m/d) 

.00384 

.00293 

.00293 

Rainfall 

a\,  (Xl0-ftm7rf') 

9.81 

1.76 

1.76 

CVf 

2.58 

1.43 

1.43 

(d) 

.6 

.6 

.6 

K*G  (m/d) 

5. 

5. 

.7 

Conductivity 

1 

1 

1 

A„  (m) 

1 

1 

1 

#000 

.1008 

.0980 

.216 

Table  4.5:  Deterministic  properties  of  the  soil  moisture  release  curve  for  the  fine  sandy 
and  sandy  loam  soils. 


Property 

Fine  sandy  soil 

sandy  loam  soil 

.35 

.496 

Or 

.063 

.11 

A 

3.7 

1.6 

K 

-.82  m 

-.85  m 

OIL 

.125m 

.125m 

0>t 

.125m 

.125m 
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i     '     i     1 i     1 i     1     r     '  I 
0.0       0.2       0.4        0.6       0.8  1.0 

Normalized  pore-water  velocity  autocovariance 


J 


Figure  4.1:  Evolution  of  the  normalized  pore- water  velocity  autocovariance  for  an 
inert  solute  in  a  heterogeneous  fine  sandy  soil  subject  to  the  Tampa  rainfall  pattern 
when  the  times  of  observation  are  adjusted  to  the  depths  of  observation. 
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Figure  4.2:  Effect  of  soil  heterogeneity  (ofnK.  =  1.  and  =  1)  in  a  fine  sandy  soil 
on  the  statistics  of  a)  the  travel  time  over  depth,  b)  the  solute  flux  over  time  at  2 
meters,  and  c)  the  solute  flux  over  time  at  6  meters. 
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Figure  4.3:  Effect  of  water  uptake  in  a  fine  sandy  soil  with  psx  =  1  on  the  statistics 
of  a)  the  travel  time  over  depth,  b)  the  solute  flux  over  time  at  2  meters,  and  c)  the 
solute  flux  over  time  at  3  meters. 
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#    Standard  deviation,  60%  water  uptake 

 Standard  deviation,  no  water  uptake 

X    Mean,  60%  water  uptake 
 Mean,  no  water  uptake 


Figure  4.4:  Effect  of  water  uptake  in  a  fine  sandy  soil  with  psx  —  —  1  on  the  statistics 
of  a)  the  travel  time  over  depth,  b)  the  solute  flux  over  time  at  2  meters,  and  c)  the 
solute  flux  over  time  at  3  meters. 
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#    Standard  deviation,  reactive  solute 

 Standard  deviation,  inert  solute 

X    Mean,  reactive  solute 
 Mean,  inert  solute 


Figure  4.5:  Effect  of  sorption  with  pr)Kd  =  1  in  a  heterogeneous  fine  sandy  soil  on  the 
statistics  of  a)  the  travel  time  over  depth,  b)  the  solute  flux  over  time  at  2  meters, 
and  c)  the  solute  flux  over  time  at  3  meters. 
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Figure  4.6:  Effect  of  sorption  with  pvKd  =  —1  in  a  heterogeneous  fine  sandy  soil  on 
the  statistics  of  a)  the  travel  time  over  depth,  b)  the  solute  flux  over  time  at  2  meters, 
and  c)  the  solute  flux  over  time  at  3  meters. 
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Figure  4.7:  Physical  description  of  the  infiltration  experiment 
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Figure  4.8:  Comparison  of  a)  the  mean,  and  b)  the  coefficient  of  variation  of  the 
nodal  solute  flux  between  the  analytical  and  the  experimental  results  in  the  middle 
of  the  four  control  planes  for  the  uniform  sandy  loam  soil  under  the  Olympia  rainfall 
scenario. 


162 


6.0e-03 


5.0e-03  - 


4.0e-03  - 


1)  3.0e-03  - 


2.0e-03  A 
1.0e-03 
O.Oe+00 


a)  Mean 


3.0e+00 


100       200  300 
Time,  in  days 

b)  Coefficient  of  variation 


400 


100       200  300 
Time,  in  days 


400 


o 
+ 

X 


Measured 
Measured 
Measured 
Measured 
—  Predicted 

  Predicted 

—  -  Predicted 
  Predicted 


at  CP4 
at  CP3 
at  CP2 
at  CP1 
at  CP4 
at  CP3 
at  CP2 
at  CP1 


Figure  4.9:  Comparison  of  a)  the  mean,  and  b)  the  coefficient  of  variation  of  the 
nodal  solute  flux  between  the  analytical  and  the  experimental  results  in  the  middle 
of  the  four  control  planes  for  the  heterogeneous  sandy  loam  soil  under  the  Olympia 
rainfall  scenario. 
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Figure  4.10:  Comparison  of  a)  the  mean,  and  c)  the  coefficient  of  variation  of  the 
nodal  solute  flux  between  the  analytical  and  the  experimental  results  in  the  middle 
of  the  four  control  planes  for  the  uniform  fine  sand  soil  under  the  Olympia  rainfall 
scenario. 
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Figure  4.11:  Comparison  of  a)  the  mean,  and  b)  the  coefficient  of  variation  of  the 
nodal  solute  flux  between  the  analytical  and  the  experimental  results  in  the  middle  of 
the  four  control  planes  for  the  heterogeneous  fine  sand  soil  under  the  Olympia  rainfall 
scenario. 
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Figure  4.12:  Comparison  of  a)  the  mean,  and  b)  the  coefficient  of  variation  of  the 
nodal  solute  flux  between  the  analytical  and  the  experimental  results  in  the  middle 
of  the  four  control  planes  for  the  uniform  fine  sand  soil  under  the  Tampa  rainfall 
scenario. 
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Figure  4.13:  Comparison  of  a)  the  mean,  and  b)  the  coefficient  of  variation  of  the 
nodal  solute  flux  between  the  analytical  and  the  experimental  results  in  the  middle  of 
the  four  control  planes  for  the  heterogeneous  fine  sand  soil  under  the  Tampa  rainfall 
scenario. 


CHAPTER  5 
CONCLUSIONS 

5.1    Summary  and  Model  Development 

This  dissertation  focuses  on  the  stochastic  analysis  of  water  flow  and  transport 
of  inert  and  reactive  solutes  in  one-dimensional  physically  and  chemically  heteroge- 
neous soils  subject  to  a  transient  flow  regime  due  to  the  variability  of  both  rainfall 
and  water  uptake.  The  primary  goal  of  this  work  was  to  develop  and  test  a  one- 
dimensional  stochastic  model  of  flow  and  transport  that  incorporates  a  temporally 
random  rainfall,  a  spatially  random  saturated  hydraulic  conductivity,  a  spatially  ran- 
dom adsorption  coefficient,  and  a  spatially  variable,  temporally  random  water  uptake. 
The  resulting  model  produced  approximate  solutions  for  the  first  and  second  moments 
of  the  soil-water  content,  Darcy's  flux,  pore-water  velocity  and  solute  flux. 

The  stochastic  model  of  unsaturated  transient  water  flow  was  developed  in 
an  Eulerian  framework.  The  one-dimensional  0-based  Richards'  equation  was  used 
as  the  governing  flow  equation  in  a  semi-infinite  domain  with  Neuman  boundary 
conditions.  The  spatial  variability  of  the  saturated  hydraulic  conductivity  and  the 
temporal  variability  of  the  water  uptake  were  described  by  stationary  random  func- 
tions in  the  governing  equation,  while  the  temporal  variability  of  the  rainfall  was 
represented  by  a  stationary  process  in  the  top  boundary  condition.  The  resulting 
stochastic  partial  differential  equation  was  solved  by  assuming  a  converging  power 
series  for  the  soil-water  content  and  by  linearizing  the  soil-moisture  diffusivity  and 
relative  hydraulic  conductivity  about  the  zero  order  soil-water  content.  The  mean  of 
the  soil-water  content  was  approximated  by  a  non-linear  equation  that  is  solved  using 
the  Newton-Ralphson  method.  The  variance  of  the  soil-water  content  was  estimated 
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at  the  second  order  in  the  Laplace  domain  before  being  numerically  inverted  using 
a  Fast  Fourier  Transform  code.  The  statistics  of  the  Darcy's  flux  and  pore-water 
velocity  were  derived  subsequently  from  the  constitutive  relationships  that  associate 
the  soil-water  content  to  Darcy's  flux  and  pore-water  velocity. 

The  stochastic  transport  model  was  established  in  a  Lagrangian  framework 
where  the  variability  of  the  solute  travel  time  is  controlled  by  the  random  properties  of 
the  rainfall,  water  uptake,  saturated  conductivity  and  adsorption  coefficient.  Because 
the  soil  was  modelled  as  a  set  of  vertical  and  independent  streamtubes  forming  a 
time-invariant  linear  system,  the  solute  flux  at  a  given  depth  is  modeled  assuming 
a  time-invariant  transfer  function  and  a  lognormal  travel  time  probability  density 
function.  The  mean  and  variance  of  the  travel  time  were  evaluated  from  the  statistics 
of  the  pore- water  velocity. 

The  main  assumption  of  the  stochastic  flow  model  was  the  linearization  of 
the  soil-moisture  diffusivity  and  relative  hydraulic  conductivity  functions  about  the 
time  and  space  invariant  zero  order  soil- water  content.  The  assumption  holds  as  long 
as  the  soil-water  content  varies  mildly  over  time  and  therefore  does  not  significantly 
exhibit  the  non-linearity  of  the  flow  equation.  However,  the  presence  of  a  water 
uptake  component  in  the  governing  equation  violates  the  assumption  of  a  time  and 
space  invariant  soil-water  content  within  the  water  uptake  depth  since  the  zero  order 
water  content  decreases  in  a  non-linear  fashion  over  depth.  In  this  particular  case, 
a  numerical  technique  is  more  appropriate  to  solve  the  resulting  stochastic  partial 
differential  equation  of  flow  with  nonconstant  coefficients. 

The  main  assumption  of  the  transport  model  is  that  solute  transport  in  the 
vadose  zone  constitutes  a  linear  system.  This  means  that  the  impulse  response  func- 
tion must  be  linear  and  time-invariant.  These  assumptions  may  not  be  met  under  a 
random  transient  flow  regime  where  the  intensity  and  occurrence  of  the  rainfall  event 
and  the  soil  initial  condition  preceding  these  events  are  temporally  variable. 


5.2    Primary  Conclusions 

The  main  findings  of  this  dissertation  are  summarized  below: 

-  The  stochastic  unsaturated  flow  model  demonstrates  i)  the  dominance  of  the 
rainfall  variability  over  that  of  the  water  uptake  and  saturated  hydraulic  conductiv- 
ity at  shallow  depths  and  ii)  the  increasing  contribution  of  the  saturated  hydraulic 
conductivity  at  deeper  depths  in  the  variability  of  the  soil-water  content,  Darcy's  flux 
and  pore-water  velocity. 

-  The  stochastic  model  of  flow  and  transport  demonstrates  that,  for  the  cases 
studied  here,  the  uncertainty  of  predicting  the  statistical  moments  of  flow  and  trans- 
port under  transient  flow  regime  is  more  controlled  by  the  soil  texture  than  by  the 
spatial  variability  of  the  saturated  hydraulic  conductivity.  The  direct  application  of 
this  finding  is  that  a  higher  priority  should  be  given  to  identify  and  delineate  precisely 
the  type  of  soils  over  a  given  study  area  rather  than  to  sample  intensively  the  spatial 
variability  within  a  soil  unit. 

-  Compared  to  numerical  Monte-Carlo  simulations,  the  stochastic  flow  model 
performs  well  as  long  as  the  linearization  of  the  soil-moisture  diffusivity  and  relative 
hydraulic  conductivity  are  accurate.  This  condition  is  encountered  when  the  temporal 
soil- water  content  variability  is  reduced  because  of  small  rainfall  variance  and/ or  high 
storage  capacity  and  low  soil-moisture  diffusivity  of  the  soil.  The  stochastic  flow  model 
predicts  flow  statistics  that  are  better  in  heterogeneous  coarse  textured  soils  when  the 
flow  is  predominantly  vertical  than  in  heterogeneous  fine  textured  soils. 

-  The  effect  of  a  spatially  variable  and  temporally  random  water  uptake  on 
the  flow  statistics  is  shown  to  i)  increase  the  variance  of  the  soil- water  content  and  to 
reduce  that  of  the  Darcy's  flux  and  pore- water  velocity  and  ii)  to  increase  the  temporal 
correlation  scale  of  the  soil-water  content,  Darcy's  flux  and  pore-water  velocity. 

-  Compared  to  numerical  Monte-Carlo  simulations,  the  stochastic  transport 
model  performs  adequately  for  all  cases  investigated.  The  quality  of  the  stochastic 
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transport  model  in  predicting  the  mean  and  the  coefficient  of  variation  of  the  solute 
flux  breakthrough  curves  improves  over  depth.  The  model  reproduces  reasonably 
well  the  mean  solute  flux  over  depth  for  all  cases  studied,  but  tends  to  underestimate 
the  standard  deviation  of  the  center  mass,  as  a  result  of  the  first  order  approxima- 
tion. Consideration  of  water  uptake  causes  a  reduction  of  the  peak  and  an  enhanced 
spreading  of  both  the  mean  and  standard  deviation  of  the  solute  flux.  The  same  effect 
is  also  found  for  the  spatially  variable  equilibrium  linear  adsorption  coefficient. 

-  The  rainfall  variability  dominates  the  solute  transport  statistics  for  inert 
solutes  but  has  a  limited  impact  on  strongly  adsorbed  solutes.  For  a  strongly  sorbed 
solute,  the  mean  rainfall  affects  the  mean  travel  time,  however  the  rainfall  variance 
is  less  important. 

-  The  Fast  Fourier  Transform  algorithm  is  an  efficient  and  fast  numerical  tech- 
nique to  solve  partial  differential  equations  with  constant  coefficients  which  provides 
good  predictions  of  the  mean  and  the  variance  of  the  flow  and  transport  processes  at 
significant  computational  savings. 

5.3    Principal  Model  Assumptions 

The  stochastic,  unsaturated  ,  transient  water  flow  model  estimates  the  mean 
and  the  uncertainty  around  that  mean  for  the  soil-water  content,  Darcy's  flux  and 
pore-water  velocity  due  to  the  variability  of  rainfall  patterns,  saturated  hydraulic 
conductivity  and  water  uptake  patterns.  Several  assumptions  have  been  made  in  the 
development  of  this  flow  model.  Following  is  a  list  of  the  primary  assumptions  and  a 
brief  statement  of  their  significance. 

-  Water  flows  in  a  one-dimensional  direction  in  a  heterogeneous  porous  medium 
where  the  air  phase  is  assumed  to  be  immobile.  Thus,  the  model  is  more  appropriate 
for  a  coarse  textured  soil  than  for  a  fine  textured  soil  where  the  flow  is  more  likely  to 
be  characterized  as  multi-dimensional.  So  far,  the  model  accounts  only  for  vertical 
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heterogeneity  of  the  soil  properties  and  does  not  consider  horizontal  variability  of  soils. 
This  limitation  can  be  relaxed  by  extending  the  current  model  to  a  three-dimensional 
flow  in  heterogeneous  soils. 

-  The  model  assumes  that  the  soil-moisture  diffusivity  function  and  the  rel- 
ative hydraulic  conductivity  functions  are  evaluated  at  a  time-and-space-invariant 
soil- water  content.  Thus,  it  does  not  account  for  any  changes  of  those  functions  over 
time,  especially  near  the  land  surface  where  the  effect  of  rainfall  is  considerable  on 
the  temporal  variability  of  the  soil-water  content  and  neither  over  space  due  to  the 
stratification  of  soils.  Accounting  for  soil  layering  would  require  i)  to  switch  from  a 
soil-water  content  based  flow  equation  to  a  pressure-head  based  flow  equation  and  ii) 
to  solve  numerically  the  governing  equations  of  the  flow  statistics  because  the  relative 
hydraulic  conductivity  function  would  be  changing  over  space. 

-  Rainfall  and  water  uptake  are  variable  over  time  but  uniform  over  space  and 
they  have  a  constant  mean  and  variance.  Thus,  the  model  does  not  account  for  either 
any  temporally  seasonal  changes  or  any  spatial  variabilities  of  the  rainfall  and  water 
uptake. 

-  Water  uptake  removes  water  from  the  soil  in  proportion  to  a  decaying  func- 
tion of  depth  multiplied  by  a  temporally  variable  sink.  Thus,  water  sink  does  not 
distinguish  evaporation  losses  from  transpiration  losses.  Also  the  model  does  not  ac- 
count for  reduced  water  uptake  under  water-limted  conditions.  For  example,  water 
uptake  is  usually  modeled  as  being  proportional  to  the  local  soil-water  content  and 
root  density. 

-  No  surface  runoff  was  considered  in  the  current  flow  model.  All  the  rainfall  is 
assumed  to  infiltrate  into  the  soil.  Therefore,  the  model  will  tend  to  overestimate  the 
volume  of  infiltrating  water  when  runoff  occurs.  Accounting  for  runoff  would  require 

i)  to  solve  the  flow  statistics  in  terms  of  pressure-head  rather  than  soil-water  content, 

ii)  to  switch  the  top  boundary  condition  from  a  Neuman  to  a  Dirichlet  type  boundary, 
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and  iii)  to  estimate  the  probability  of  runoff  occurence  for  a  given  combination  of  soil 
and  rainfall  pattern. 

-  The  flow  model  accounts  for  the  spatial  variability  of  the  saturated  hydraulic 
conductivity  while  maintaining  the  pore  size  distribution,  the  porosity,  and  the  bub- 
bling pressure  uniformly  constant  over  space.  Because  of  the  Hagen-Poiseuille  re- 
lationship, the  most  variable  parameters  such  as  the  pore  size  distribution  and  the 
bubbling  pressure  should  be  adjusted  to  the  variations  of  the  saturated  hydraulic 
conductivity.  However,  the  spatial  variabilities  of  the  bubbling  pressure  and  pore  size 
distribution  have  little  effect  on  flow  and  transport  statistics  (Sassner,  1995,  p.  66) 
and  therefore  can  be  considered  to  be  uniform  over  space.  Note  that  varying  the 
bubbling  pressure  has  a  greater  impact  on  flow  ensemble  moments  than  a  change  of 
pore  size  distribution. 

-  The  stochastic  flow  model  considers  the  spatial  variability  of  the  saturated 
hydraulic  conductivity  with  a  maximum  variance  of  the  log-saturated  conductivity 
equivalent  to  the  unity.  While  this  value  satisfies  the  solution  technique  in  the  con- 
verging power  series,  it  is  also  representative  of  field  measurements  of  saturated  hy- 
draulic conductivity  (Gelhar,  1993,  p.  291). 

The  stochastic  advective-dispersive  transport  model  estimates  the  mean  and 
the  uncertainty  around  that  mean  of  the  solute  flux  concentration  for  inert  and  re- 
active solute  due  to  the  variability  of  rainfall  and  water  uptake  patterns,  saturated 
hydraulic  conductivity,  and  linear  adsorption  equilibrium  coefficient.  Below  are  listed 
the  assumptions  of  the  transport  model  and  their  significance. 

-  The  soil  is  represented  as  a  vertical  streamtube.  The  model  considers  only 
vertical  dispersion  along  a  streamtube.  Neglecting  the  horizontal  dispersion  means 
that  the  amount  of  solute  passing  a  given  depth  will  be  overestimated. 

-  The  adsorption  process  is  described  by  a  linear,  reversible,  equilibrium  model. 
If  sorption  is  described  by  a  non-linear  isotherm,  the  amount  of  solute  passing  a  given 
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control  plane  will  depend  upon  the  solute  flux  concentration.  If  adsorption  is  not 
instantaneous,  the  model  will  tend  to  underestimate  the  amount  of  chemical  passing 
a  given  depth.  If  adsorption  is  not  irreversible,  the  amount  of  chemical  passing  a 
given  depth  will  be  overestimated. 

-  The  half-life  time  that  measures  the  biodegradability  of  the  chemical  is  not 
considered.  Thus,  the  model  will  overestimate  the  amount  of  chemical  passing  a  given 
depth.  However,  the  half-life  time  can  be  integrated  in  the  impulse  response  function. 

-  The  impulse  response  function  is  assumed  to  be  time-invariant  and  calibrated 
for  a  given  soil- water  content  profile  in  the  soil.  Thus,  it  does  not  account  for  the 
temporal  variability  of  the  soil-water  content  profile  at  the  time  of  solute  application. 
Depending  upon  that  profile,  the  transport  model  will  either  underpredict  the  solute 
mass  passing  a  given  depth  if  the  soil  was  initially  nearly  saturated  or  overpredict  the 
solute  mass  passing  a  given  depth  if  the  soil  was  initially  dry. 

5.4    Applications  and  Further  Research 

The  soil-water  content  statistics  developed  in  this  dissertation  should  be  use- 
ful for  the  parameterization  of  land-surface  models  in  order  to  improve  predictions 
of  land-surface  water  flux  interactions  and  global  climatic  changes  (Budyko,  1956; 
Delworth  and  Manabe,  1988). 

The  transport  model  can  be  used  in  association  with  a  Geographical  Infor- 
mation Systems  (GIS)  in  order  to  provide  field-scale  risk  assessment  of  nonpoint 
source  pollutant  in  the  vadose  zone.  The  three-dimensional  soil  would  be  represented 
by  a  cluster  of  vertical  and  independent  soil  columns  over  which  the  current  one- 
dimensional,  stochastic,  advective-dispersive,  transport  model  can  be  applied.  No 
exchange  of  solute  mass  will  occur  between  those  soil  columns.  A  recent  paper  by 
Corwin  et  al.  (1997)  discusses  the  possibility  of  coupling  flow  and  transport  mod- 
eling to  GIS  in  order  to  model  nonpoint  source  pollutant  in  the  vadose  zone.  This 
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approach  is  emphasized  by  Jury  (1996)  who  suggests  that  a  stochastic  stream  tube 
model  seems  the  most  compatible  with  GIS  because  "it  utilizes  a  relatively  simple 
local  process  driven  by  parameters  that  might  be  associated  with  soil  morphological 
features,  and  could  be  integrated  up  to  a  large  scale  by  simple  arithmetic  averaging 
over  the  local  sites" .  The  current  one-dimensional  stochastic  transport  model  is  well 
adapted  for  this  objective  that  has  already  been  met  for  the  deterministic  pesticide 
fate  model,  CMLS  (Nofziger  et  ai,  1987),  by  Foussereau  et  al.  (1993). 

The  current  stochastic  transport  model  should  also  be  useful  to  predict  solute 
flux  statistics  in  the  coupled  soil-groundwater  system  described  by  Destouni  and 
Graham,  (1995).  In  this  approach,  the  solute  transport  in  both  zones  is  modeled  by  a 
response  impulse  function  multiplied  by  a  input  function  and  the  outputs  coming  from 
the  vadose  zone  correspond  to  the  inputs  of  the  saturated  zone.  Therefore,  following 
this  methodology,  the  outputs  of  the  current  stochastic  transport  model  can  be  used 
as  the  input  function  of  the  transport  model  in  the  saturated  zone  in  order  to  predict 
solute  flux  statistics  in  groundwater. 

The  capabilities  of  the  current  model  could  be  improved  by  incorporating 
multi-dimensional  flow  and  also  by  incorporating  a  water  uptake  flux  as  a  function  of 
the  local  soil-water  content.  The  rainfall  and  the  water  uptake  series  were  assumed  to 
be  second  order  stationary.  In  order  to  be  more  realistic,  these  restrictions  could  be 
relaxed  by  considering  stationary  variances  around  temporally  variable  means  that 
follow  seasonal  patterns. 
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APPENDIX 
LAPLACE  TRANSFORM 

A  detailed  and  concise  review  of  the  relationship  between  the  Fourier  and 
Laplace  transforms  is  presented  by  Brigham,  (1988,  p.  200).  Here  are  summarized 
the  major  points.  The  Laplace  transform  of  a  real  function  of  time,  g(t),  is  given  by 

r+oo 

G(si)=  /      g(t)e-^dt  (A.l) 
Jo 

If  Si  =  C\  +  i2irfi  then  the  previous  equation  becomes, 

r+oo 

G{c1+i2nf1)=  /      [g(t)e^]e-l2^dt  (A.2) 
Jo 

where  c\  is  a  real  constant,  /i  is  the  frequency  and  i  is  the  complex  number 
equivalent  to  If  g(t)=0  for  t  <  0  then  the  lower  limit  of  equation  (A.2)  can  be 

replaced  by  — oo.  Hence,  the  previous  equation  can  be  recast  in  term  of  the  Fourier 
transform, 

^)e-ci*<^G(c1  +  z27r/1)  (A.3) 

Therefore  the  inversion  of  the  Laplace  transform  consists  of  i)  rewriting  G(si) 
by  replacing  the  variable  sx  by  cx  +  i2nfi,  ii)  inverting  G{cx  +  i2nfi)  by  using  a  Fast 
Fourier  Transform  algorithm  and  iii)  multiplying  the  inversion  result,  g(t)e~Clt,  by  eClt 
in  order  to  recover  g(t).  The  same  procedure  is  applied  for  2-D  Laplace  transform, 

g(t,  Oe-Clf-C2t'  ^  G{Cl  +  i2irf\ ,  c2  +  i2nh)  (A.4) 
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